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NUMERICAL  SOLUTION  FOR  THE  SUPERSONIC  INVISCID  FLOW  AROUND 
NON-SYMMETRICAL  BODIES 


by  Guo  Zhiquan  and  Zhang  Lumin 

(Aerodynamic  Research  and  Development  Center  of  China) 


Abstract 


This  paper  describes  a  numerical  method  for  evaluation  of 
the  three-dimensional  inviscid  supersonic  flow  around  the  re¬ 
entry  vehicles  of  non-symmetrical  bodies.  The  method  is  devoted 
to  determining  the  initial -boundary  value  problem  for  three 
independent  variable  first  order  stationary  quasi-linear  hyper¬ 
bolic  partial  differential  equations.  MacCormack 's  explicit 
finite  difference  schemes  with  second  order  accuracy  is  adopted 

The  numerical  examples  of  the  bent  sphere  cones  and  the  re 
entry  vehicles  with  control  flaps  have  been  worked  out  to  demon 
strate  the  capability  of  this  computation  method.  The  accuracy 
of  this  technique  has  been  found  to  be  satisfactory. 


Symbols 

A 


a 


a 


e' 


b 

e 


C 

P 


J 

o 


K 


the  total  drag  on  the  body's  surface  in  front 
of  the  solution  section 

local  sonic  velocity  (  a  =  VfFT?) 

the  curvilinear  -  ellipsoid  long  and  short  semi- 
axes  of  the  body's  surface  on  the  solution 
section  (see  figs.  1  and  2) 

pressure  coefficient  on  surface  (c  =  P  —  I 

'  1 r  *1 

total  number  of  radial  and  peripheral  solution 
crunodes  in  the  solution  section 


r,z,V 


r.  ,r 

D  S 


rb  'rs 
z  : 


local  and  free  flow  Mach  numbers 
pressure 

radial,  axial  and  peripheral  (meridian  angle) 
coordinates  in  the  cylindrical  coordinate  system 

the  heights  of  the  body's  surface  and  shock  wave 
surface  (radial) 

s<p  the  first  order  derivatives  of  or  rg  for 
the  lower  angular  coordinates 


3^  the  second  order  derivatives  of 
”  lower  angular  coordinates 

flow  around  blunt  nose  radiup  of  body 


for  the 


S  entropy  (S=P /f>Y) 

u,v,w  axial,  radial  and  peripheral  velocity  components 

u^  ,  v^  (W^  axial,  radial  and  peripheral  components  of  free 
flow  velocity 

Vneo  normal  velocity  in  front  of  shock  wave 

a  free  flow  attack  angle 

y  specific  heat  ratio,"/ =1.4  in  ideal  gas 

&Z.  axial  radial  and  peripheral  network  step  lengths 

$  included  angle  (below  called  bent  angle)  of  front 

N  and  back  cone  axes 

£  included  angle  of  the  tapered  plane  and  body 

axis  (called  "tapered  surface  angle"  for  short) 
in  the  "tapered  cone" 

£ _  included  angle  (control  angle)  of  the  control 

wing  and  body  axis 

&C'  ^FC  semi-cone  angle  of  each  cone 

i  (h  radial  and  peripheral  coordinates  after  trans- 

>  *  formation 

K'  first  order  derivatives  of  §  for  lower  angular 

r  coordinates 

P  density 

derivative  of  #  for  ? 


Upper  Corner  Symbols 

n  (axial)  solution  section  code 
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Lower  Corner  Symbols 

i  code  of  sectional  cylindrical  coordinate  system 

j,k  radial  and  peripheral  codes  of  solution  crunodes 

I.  Preface 

The  great  majority  of  high  performance  maneuverable  re¬ 
entry  vehicles  draw  support  from  the  non-symmetricality  of  the 
configuration  to  attain  control  moment  and  to  cause  it  to 
produce  a  predetermined  balanced  angle  of  attack.  In  this  way, 
there  is  then  maneuvering  flight.  Therefore,  attention  must  be 
given  to  the  solution  of  the  three-dimensional  irtviscid  flow 
around  non-symmetrical  bodies. 

The  flow  field  of  the  supersonic  inviscid  flow  can  be 
divided  into  two  sections  according  to  differences  in  quality: 
the  subsonic  and  transonic  area,  and  the  rear  supersonic  area. 
This  paper  only  discusses  the  solution  of  the  rear  supersonic 
area  and  the  nose  flow  fields  are  taken  from  references  [3,4]. 

Solutions  for  the  rear  area  of  non-symmetrical  bodies  began 
during  the  1960's  in  foreign  nations  and  during  the  1977-1979 
period  many  technical  achievements  were  attained.  By  the  end  of 
1979,  domestic  work  showed  results  as  seen  in  references  [1,2] . 
The  major  distinctions  of  this  paper  are: 

(1)  The  coordinate  systems  are  different.  For  bent  sphere 
cones,  different  oriented  cylindrical  coordinate  systems  are 
used  in  different  cone  regions  yet  they  are  not  the  coordinate 
systems  where  "the  original  points  of  the  coordinates  shift 
along  the  body's  axis  and  the  direction  is  invariant." 

(2)  The  beginning  equations  are  a  set  of  conservation  law 
flow  equations  whose  differential  forms  are  separately  given  in 
reference  [2] . 


(3)  There  are  differences  in  the  solution  method  for  the 
body 1 s  surface . 

(4)  When  there  was  a  relative  lack  of  experimental  and 
numerical  computations,  we  used  the  total  mass,  total  energy 

and  total  momentum  conservation  examination  technique. 

\ 

t 

The  method  used  in  this  paper  can  be  extended  to  situations 
where  there  is  non-symmetrical  yaw  (slip  angle)  and  the  results 
obtained  were  satisfactory  (to  be  published  in  a  separate 
article) . 

We  would  like  to  thank  Zhang  Hanxin,  Li  Songbo  and 
Lao  Shuchun  for  all  of  their  help  in  the  work  for  this  paper. 

II.  Basic  Equations  and  Boundary  Conditions 

1.  The  Selection  of  the  Coordinate  Systems1 

The  bent  sphere  cones  (including  the  multiple  region  bent 
vehicle)  were  mostly  multiple  region  cones  and  the  body's  axis 
was  not  on  the  same  straight  line.  Therefore,  different  body 
axis  cylindrical  coordinate  systems  oi~ziri  9*^  were  used  for  the 

different  cone  regions.  Given  this  type  of  regional  cylindrical 
coordinate  system,  each  region's  corresponding  basic  equation 
had  completely  identical  forms.  To  be  sure,  the  joining  areas  of 
the  front  and  back  coordinate  systems  could  present  trouble: 
firstly,  as  the  beginning  surface  Zi+1=Zi+1  0  °f  the  back  axis 

coordinate  system  computation,  we  must  give  the  "initial  value"; 
secondly,  as  the  end  region  (Zg^Q  ^  Z^  ^  Zeil^  tlie  ^ront  axis 


1.  Zhang  Hanxin  offered  valuable  suggestions  for  the  determination 
of  the  coordinate  systems. 
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coordinate  system  computation,  it  not  only  has  aerodynamic 
solution  integral  limitation  problems  but  it  also  only  provides 
interpolation  data  results  for  the  above  mentioned  "initial 
values."  These  are  easily  resolved. 


w 


Fig.  1(a)  Coordinate  System  and  Sketch  of  a  Bent  Sphere  Cone 

Key:  1.  Beginning  section  I 

2.  Beginning  section  II 

3 .  Shock  wave 

4 .  Shock  wave 


( b ) 

Fig.  1(b)  Contour  in  Sections  A-A1  and  B-B' 

2.  Basic  Equations 

In  the  (body’s  axis)  cylindrical  coordinate  o-zrV  (omitting 
the  lower  angular  coordinate  i) ,  if  the  body’s  surface  and 


shock  wave  surface  are  separately  indicated  as 


r  ”r*(  2  ,  CP),  r«r,<2,  <p  ) 


(1) 


and  we  introduce  the  following  transformations 


2  =  2, 


<t>=tg"(  h;t 


(2) 


(for  ag  and  see  fig.  1(b)),  using  this  type  of  non-dimension¬ 
al  izat ion:  P  and  fi  are  separately  divided  by  free  flow  pressure 
P»  and  free  flow  density  P^  ;  the  velocity  is  divided  by 

V  P m/P*  z  ;  r  etc.  indicates  that  the  length  is  divided  by 


blunt  nose  radius  R^.  Thus,  the  three-dimensional  inviscid 
stationary  flow  equations  are 


*A  . 

*2  di 


*+»l+wl+~' =C, 


(3) 

(4) 


In  the  equations 


C,=  YM1+— 2J1 

X=*(P«,  Pul+p,  Puv,  Puts))1 

3-4.A+4rB*+4»c* 


and 


Pu 

•Pie  ' 

Po 

B*- 

Puv 

Pv*+p 

C*«— i- 
,  c  — 

Puie 

Pulo 

•  5*-r~ 

Puo 

P(l»*  — ifl*) 

,  Pvw 

.Puj*+  P  - 

.  2Pute  < 

The  region  of  solution  is 

0  1  *  0  ji  (  i  ■  1|  2  y  **•) 

/ 

3 .  Boundary  Conditions 

(1)  On  the  body's  surface  (  ^»  =0),  there  are  non -permeating 
conditions : 


urt>—  v  +uT»,/r»—  0 


C5) 


(2)  On  the  shock  wave  surface  (£»  =1),  there  is  the  Rankine 
Hugoniot  relational  formula.  It  can  be  rewritten  as 


v 

-  2*2 

p  -vi./ ( i  +r:.-p) 

r'.  “  — —  <«-(*-  -  w-o/r,) 

+K„i/(w.-«;.r,yr0,4  (ui-f'l.X  1 

«  +  - i-)/  y/  i  -t-r'Jr2, 

«  =u.-(u-o.)r,i 
»=».-(»—  Vm)r,Jr, 


(6) 
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s*.  •-  ««i»  «*? 


In  the  equations,  Vn#  is  the  normal  velocity  in  front  of  the 
wave  and  .v#  and  vm  are  separately  the  (free  flow) velocity 
components  in  front  of  the  wave : 


«-  =■  ^  Y  M-coa  o 

vm=*  —  V  y  \j_sin  a cos<t>  ’ 

wm  sm*'  Y  M.sin  a  sin  <J>  J 


(7) 


(3)  If  we  take  into  consideration  the  problem  of  maneuver- 

/ 

ing  within  the  surface  and  that  the  body  has  a  symmetrical  sur¬ 
face  (  $  =0,7 r  )  we  can  apply  symmetrical  conditions  on  it. 

4.  Initial  Conditions 

On  the  left  end  surface  (initial  section) z . =z .  of  the 

1  10 

solution  region,  the  initial  values  of  its  flow  field  and  shock 
wave  parameters  car.  be  interpolated  from  the  values  of  the 
upstream  area. 

For  the  initial  boundary  value  problem  of  the  set  of  hyper¬ 
bolic  equations  (3) ,  its  applicable  conditions  are  the  initial 
value  sections  of  the  space  towards  the  curved  surface. 

III.  Differential  Equations  and  Boundary  Processing 

If  the  equidistance  on  the  variation  region  [0,1]  of 
gives  crunodes,  the  distance  between  the  crunodes  is 
A^o  =(Jo-l)-S  if  the  equidistance  on  the  variation  region  of 

gives  Kq  crunodes,  A$  =  'Tf/i  Kq-1) .  In  the  solution  of  each 

layer  along  the  axis,  after  the  network  point  flow  parameters  on 
the  z=z{  }  section  are  known,  all  of  the  points  on  the 

z=z  ^n+1  ^  =z  ^  +  A  z  section  can  be  separately  calculated  as  follows 
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1 .  Calculation  of  Inner  Points 


In  the  network's  inner  points  (l  ^  j  K.  Jq,  1  K  k  <  Kq)  , 

MacCormack's  explicit  finite  difference  schemes  with  second  order 
form  are  used  for  the  structure  of  the  set  of  equations  (3) : 


Am  -  a*.*-  b7.*)/AI c;,»)/AO + C  8 ) 

"ti “  — ^ *+■  ^4 /»* ~  B  !>*—  B  ;-i>*)  /A4 

+(c;,*~  C;.*-i)/ A<t>  +  Z);,*]}  (8a) 


Equation  (8)  is  the  predicted  step,  is  the  predicted  value 

and  equation  (8a)  is  called  the  correction  step. 

After  each  step  of  computation  attained  A  .  the  simul- 

J  .  * 

taneous  Bernoulli  integral  can  then  obtain  the  predicted  (or 
corrected)  values  of  the  flow  field's  inner  point  flow  parameters 


2 .  Calculation  of  Shock  Wave  Surface 

At  this  time,  it  is  also  necessary  to  determine  shock  wave 

parameters  r  ,  r  and  r  .  For  this  reason,  the  set  of 
s 

z 

equations  in  (6)  are  not  closed.  We  can  simply  use  a  method 
analagous  to  calculating  the  inner  points  (yet  in  the  predicting 
step,  we  change  to  use  the  backwards  difference  quotient  for  the 
differential  quotient  of  £  ) .  First  we  determine  the  wave  and 
then  the  prediction  or  correction  value  of  pressure  PjQ  As 
for  rs,  we  use  the  improved  Youla  (?)  broken  line  method;- 
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Predicted  step 


(9) 


- (r  Jv=6-  )  J+  4  z(r  )  JJ 


s'k  'rs  •  \C 
n+1 


V  k 


Correction  step  —  (r  ^  k  =  §  [<r8^  *(?,)  Jilr,  1„] 

z 

fr*  is  determined  by  using  a  central  difference  quotient 


(9a) 


(rO*-|_i__CCr }  3 


k  ~  l ,  K, 

1  <  A  </f. 


(10) 


In  this  way,  we  can  compute  the  flow  parameters  after  the  wave 
from  equation  (6)  and  the  prediction  values  and  correction  values 
of  the  shock  waves - 

3 .  Computations  of  Body ' s  Surface 

On  the  body's  surface,  aside  from  the  boundary  conditions 
in  equation  (5)  the  flow  should  satisfy  the  following  quasi-wave 
characteristic  equation 


a  (G  -arti) 

l 

frh  (no*- A  rf-g1 

1  aP 

dz 

1 1:  — a* 

I4,  *4 

r»<; 

l  ’  \  o  G  —  art'  rk  ) 

'  d»t> 

,T--*r  -Tr)- *•{-£-  .-SO 

+  Pm»jj  4-  -SJ.  [(«■„_- anjn+m+jr.)] 


(ID 


Cl  1(« 1 * G  «  ✓?(  1  +rl'+rl~jrl)-ai(  1  +  ri,/rj)) 


Key:  1.  In  the  equation 
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and  flow  characteristic  equations 


_ »_ 

6z  r, u  6<t> 


(12) 

(13) 


We  also  used  MacCormack ' s  explicit  finite  difference 
schemes  with  second  order  form  to  calculate  equations  (11) -(13). 
When  we  considered  the  body's  surface  j  =  l,  the  differential 
quotients  of  £  were  all  replaced  by  difference  quotients  when 
in  the  first  and  second  order. 

It  must  be  pointed  out  that  although  the  above  method  of 
computing  the  body's  surface  was  able  to  accurately  reflect  the 
inviscid  flow  characteristics,  yet  for  the  computation  of 
longer  cone  sections  (for  example,  when  its  length  is  greater 
than  50  times  the  nose  radius)  and  the  leeward  surface  it  cannot 
go  down  continuously.  This  is  because  as  the  rear  body  high 
entropy  layer  becomes  thinner,  the  local  Mach  number  on  the 
body's  surface  becomes  smaller.  As  soon  as  the  axial  velocity  is 
subsonic,  the  applicable  conditions  of  the  initial  boundary 
value  problem  of  the  hyperbolic  equations  are  destroyed.  For  the 
the  computation  to  go  down  continuously,  the  entropy  must  be 
decreased.  We  can  begin  from  any  suitable  section  and  change  and 
use  the  set  of  equations  in  (3) ,  (4)  and  (5)  (deleting  the  third 
equation  in  the  equations  of  (3)).  The  large  majority  of  exper¬ 
iment  computations  have  shown  that  this  type  of  non-strict 
decreased  entropy  method  was  functional  and  the  computation 
results  of  the  body's  surface  pressure  were  quite  accurate. 

4.  Handling  the  "Singular  Line"  of  the  Body's  Surface 

On  the  "edge"  of  the  non -continuous  slope  of  the  body's 

11 


surface,  the  flow  is  multi-valued.  The  non-continuity  of  the  flow 
only  occurs  on  the  normal  plane  of  the  edge  line.  The  points  on 
the  edge  are  "expansion  angle  points"  or  "compressed  angle 
points."  They  are  determined  by  whether  the  air  flow  deflection 
angle  on  the  edge  line  normal  plane  of  the  point  is  positive  or 
negative  (if  the  body's  surface  curve  on  the  normal  plane  is 
protruding,  then  the  air  flow  deflection  angle  is  defined  as 
"positive",  otherwise  it  is  "negative").  After  calculating  the 
upstream  flow  value  of  a  certain  computed  point  on  the  edge,  we 
used  the  method  found  in  reference  [  1]  to  determine  the  down¬ 
stream  flow  value  of  the  point.  However,  there  is  an  error  in 
equation  (4.9)  in  reference  [1]  and  it  should  be  changed  to 


X*-f  c,X,+etX4-c,*«  0  1 

c,--(  1  +2/M‘)-Ysin,e1 

e,» (2M*  +  l)/M‘-+-((  Y  +  l  )*/4  +  (  Y  —  1  )/M*)»n*0,  (14) 

c,- -cosfy/M4 

Here,  X  is  the  sine  square  of  the  shock  wave  angle  and  Q i  is 
the  air  flow  deflection  angle. 

IV.  Stability  Conditions 

Step  length  z  must  be  limited  by  the  stability  conditions. 
A  necessary  condition  for  MacCormack's  explicit  finite  difference 
equations  is  that  for  any  one  point  in  the  solution  region,  the 
dependence  region  of  the  differential  equation  should  be  com¬ 
pletely  included  in  the  dependence  region  of  the  difference 
equation,  that  is 

max  ((°i );,»,  (ai )/,*»  (CT» )m)  05) 

Kk<K, 
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In  the  formulas 


o,  -  -  ^",1- jj  *  (4. + uir  +  utU/r)  -  a*U 

+  a  v/(ot,4-u»4*/r)’  +  (u‘ — a1) (4J  +  WM)  j 

°* — rvfcsr 

o, -jrry jl  <  [»6.+»«,+-7-{5.+*.-||-)] -«,U 

+ •  /K+^.+t.-f!)]' +<“,-“,)[5!  +(6-+*-i»  J/rl  11 

Because  of  this  ' 

A-z'^Afc/max  <(ai )>,*t  (°»)/.*}  (16) 

x<  i  <Js 

l<k<K, 

JuU  can  generally  take  0.9 


V.  Computation  Examples 

1.  cone  -  Cylinder  -  Flare  Symmetrical  Bodies 

We  programmed  and  computed  the  bent  sphere  cone  (€^=0) .  It 

can  be  seen  from  fig.  2  that  the  two  types  of  computations  are 
basically  the  same  and  match  well  the  experimental  results. 
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Fig.  2  0^  Curves  as  Compared  to  Computational  Results  of  [5] 

and  the  Experimental  Results  of  [6] 

Key:  1.  Experimental  results  of  [6] 

2.  Confutations  of  this  paper 

3.  Computations  of  reference  [5] 

2 .  Bent  Sphere  Cone 

The  flow  fields  were  computed  for  the  configurations  shown 
in  fig.  1(a),  6FC=10°,  @B=5°,  under  different  bend  angles 

(  £  N= 0°  ~  6°),  different  free  flow  Mach  numbers  (M^  =5,8,10) 

and  different  attack  angles  (a=0°~  10°).  Fig.  3  gives  the  partial 
pressure  distribution  curve  and  fig.  4  presents  the  pressure  of 
the  different  bend  angles  )  which  change  the  curve  in 

accordance  with  the  attack  angle  and  have  similar  variation 
tendencies  as  reference  [2]  (in  this  reference,  they  only  give 
the  curve  but  no  measurements  so  that  there  is  no  way  of  making 
a  quantitative  comparison) . 
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Fig.  4  /L  Vs.  a  Curves  of  a  Bent  Sphere  Cone 

«»  "8.  0FC'1<>O.  V5°-ZU=9-6'ZL2*20 

3 .  Tapered  Cone  With  Control  Flaps 
Its  configuration  is  shown  in  fig.  5. 
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Fig.  5  Sketch  of  a  Tapeted  Cone  With  Control  Flaps 

We  computed  different  free  flow  Mach  numbers  Ma?  =5  and  10, 
different  attack  angles  <x=0°  10°  and  different  body  shapes 

(tapered  on  a  single  surface,  on  two  surfaces,  on  four  surfaces 
etc.)  for  @c=10°,  4^=7. 5°  and  ?>2®14.5  ^21°.  Fig.  6  gives  a 

set  of  pressure  distribution  curves. 


* tl  ‘LI  * 


Fig.  6  Pressure  Distribution  of  a  Tapered  (on  Both  Sides) 
Cone  With  Control  Flaps 

4.  Tapered  Double  Cone 

This  is  the  bent-tapered  cone  (  $  N=0)  and  the  only- 
complete  results  found  in  the  references.  This  paper  uses  the 
"bent-tapered  cone"  sequence  to  carry  out  computations  and  a 
comparison  with  the  results  of  reference  [1]  are  derived  from 
fig.  7  and  listed  in  the  following  table: 
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■ 

S  i 
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xstcn 
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0.3569 

0.2576 

0.3357 

0.3362 

0.2238 

0.2246 

IS  J 

0.4402 

0.4399 

0.4079 

0.4074 

0.2336 

0.2331 

j4.i  ; 

0.0746 

0.0744 

0.0646 

0.0643 

0.0196 

0.0195 

3S.6  ! 

0.0497 

0.0354* 

0.0645 

0.0642 

0.0195 

0.0194 

31.6  | 

0.0347 

0.0363* 

0.0625 

0.0624 

0.0193 

0.0193 

40.6  1 

0.0362 

0.0381* 

0.0599 

0.0603 

0.0193 

0.0193 

43.1  j 

0.0375 

0.0378 

0  0561 

0.0563 

0.0194 

0.0194 

44.6 

0.0373 

0.0375 

0.0536 

0.0536 

0.0193 

0.0194 

(.12?*'”  ’30’’ 


(7)h-  ft  *  m 

ft  X 
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At  *  13 

/o  *  21 

K,m  IS 

/.-is  4 

K,-  25 

/o“21 

Ko  ■  15/  q  \ 

(ftB) 

/••IS 

K.-3S 

C0*,96*Drt25t 

C9S\  *3(11 

Table  1 

Key:  1.  This  paper 

2.  Reference  [1] 

3 .  This  paper 

4.  Reference  [1] 

5 .  This  paper 

6.  Reference  [1] 

7 .  Computation  network 

8 .  This  paper 

9.  Equidistance 

10.  Reference  [1] 

11.  25  equidistant  points  in  10,96°],  8  non- 
equidistant  points  in  [96°,  "If] 

12.  The  values  of  SP  =30°  and  90°  in  this 
table  [ 1]  are  interpolations  made  by 
the  authors 


There  is  only  a  0.2  percent  difference  between  the  two  in 
front  of  the  "tapered  surface."  The  difference  within  the  "tapered 
surface"  is  somewhat  larger  and  the  largest  difference  occurs  on 
the  extreme  left  end  of  the  "tapered  surface."  The  reason  for 
this  is  that  reference  ( 1]  does  not  use  equal  entropy  relational 
formula  (13)  in  computing  the  body's  surface  but  rather  uses  the 
so-called  decreased  entropy  method  (the  entropy  attained  by 
interpolating  the  inner  point  entropy  on  the  body's  surface).  In 
two  types  of  computations,  the  upstream  pressure  value  on  the 
left  end  point  of  the  "tapered  surface"  (z=35,  V  =0°)  were 
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identical.  Due  to  the  decreased  entropy,  the  point's  upstream 
density  was  higher  than  the  value  of  the  equal  entropy  method 
(reference  [  1]  does  not  give  density  values  and  thus  no  compar¬ 
ison  can  be  made) ,  and  the  local  Mach  number  was  higher  and 
because  of  this  expanded  sharply.  Therefore,  the  point's  down¬ 
stream  pressure  value  was  also  lower.  The  differences  in  the 
initial  end  of  the  "tapered  surface"  which  was  caused  by  the 
different  handling  of  the  entropy  had  a  very  small  effect  on  the 
pressure  of  the  body's  surface  outside  the  "tapered  surface." 
Even  inside  the  "tapered  surface"  its  differences  also  gradually 
decreased  and  disappeared  in  accordance  with  the  downstream 
propulsion . 


Fig.  7 
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Fig.  7  Pressure  Distribution  of  a  Tapered  Double  Cone  as 
Compared  With  the  Results  of  Reference  [  1] 

Key:  1.  Computations  of  this  paper 

2.  Computations  of  reference  [1] 

VI.  Analysis  of  Precision 

Given  that  there  is  a  lack  of  comparison  of  experimental 
and  other  numerical  computations,  to  examine  the  reliability  of 
this  method,  we  mainly  used  the  total  mass,  total  energy  and 
total  momentum  conservation  technique  to  carry  out  investiga¬ 
tions  .  ^ 

The  shock  wave  on  the  solution  section  was  taken  as  the 
outward  boundary  to  make  a  flow  tube.  If  the  outer  flow  control 
surface  is  s  ,  the  annular  region  between  the  shock  wave  and 

body's  surface  on  the  solution  surface  is  The  total  mass  and 

total  momentum  are  recorded  as  Gq  and  Mq,  and  the  corresponding 

symbols  on  Sj  are  G1  and  m1 .  Furthermore,  if  the  body  surface's 

total  drag  in  front  of  the  solution  section  is  A,  then  the  total 
mass  and  total  momentum  conservation  can  be  written  as: 


Mt  =  Mt  +  A 


(17) 

(18) 


1.  Zhang  Hanxin  and  Gao  Shuchun  participated  in  discussions  of 
this  problem. 
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With  Control  Flaps 


Fig.  9  Errors  of  Computed  Total  Momentum  of  a  Bent  Sphere 
Cone 

Computations  show  that  for  a  scattered  network,  the  pre¬ 
cision  is  lower  (e.g. ,  Jo=ll,  the  max (4  M)  of  the  bent  sphere 

cone  is  about  1  percent,  the  tapered  cone  is  about  0.5  percent); 
the  precision  gradually  rises  as  the  network  becomes  denser  and 
when  Jq=31  the  bent  sphere  cone  is  0.1  percent  and  the  tapered 
cone  is  about  0.3  percent. 

As  shown  in  figures  8  and  9,  the  larger  errors  occur  in  the 
initial  section  and  extreme  left  end  of  the  "tapered  surface". 
This  is  due  to  the  use  of  the  interpolation  formula  in  these 
sections  which  causes  larger  truncation  errors. 

To  sum  up,  this  method  can  attain  relatively  high  Drecision 
under  total  mass,  total  energy  and  total  momentum  conservation. 
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EXPERIMENTAL  STUDY  ON  HEAT  TRANSFER  TO  ROUGH  WALLS 


by  Tian  Wenbing,  Li  Xue  and  Mao  Mingfang 
(Beijing  Institute  of  Aerodynamics 

Abstract 

Experimental  results  of  heat  transfer  to  rough  walls  are 
given  in  this  paper  for  sphere  cone  models  at  Mach  number  5. 

The  nose  radius  of  the  models  is  27.4  mm  and  its  base  diameter 
is  60  mm.  Five  models  have  been  tested  and  the  different  rough¬ 
ness  in  its  bend  dianeters  are  separately  0,0. 3, 0.5, 0.7  and  0.9 
mm.  The  tests  were  conducted  in  a  conventional  hypersonic  wind 

tunnel  at  total  pressures  p  of  10-45  kg/cm  and  Reynold’s 

^  £ 

nunbers  Rep  of  (0.8-3. 6)  x  10  • 

The  test  results  indicate  that  the  smooth  wall  model  heating 
is  the  laminar  flow  heating,  its  heat  flux  stagnation  point  is 
quite  close  to  the  laminar  flow  theoretical  data,  and  the 
influence  of  roughness  at  low  total  pressure  (lOkg/cm  )  occurs 
mainly  to  promote  the  transition  and  development  of  a  boundary 
layer.  With  the  increasing  total  pressure  in  the  wind  tunnel  the 
static  pressure  on  model  and  local  Reynolds  numbers  correspond¬ 
ingly  increase.  In  this  case,  the  effect  of  roughness  on  heat 
transfer  becomes  remarkable  and  the  most  remarkable  region  appears 
at  the  sonic  point  region  on  the  nose.  At  the  highest  total 
pressure  (pfc=45  kg/cm  )  and  with  the  largest  roughness  (d=0.9  mm), 

the  ratio  of  rough  wall  heat  flux  to  laminar  flow  smooth  one 
(qi/qyo)  could  be  up  to  4  except  at  stagnation  point  where  it 

could  approach  6.  Its  raise  seems  to  be  related  to  the  local  shape 
change  in  the  vicinity  of  the  stagnation  point. 
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I .  Preface 


The  problem  of  roughness  and  equivalent  roughness  was  already 
studied  during  the  1930's  [1].  Recently,  following  the  extended 
use  of  the  warhead  ablation  plan,  this  has  been  given  greater 
attention.  During  the  re-entry  process,  the  sonic  point  region 
on  the  nose  becomes  very  rough.  For  example,  the  average  height 
of  these  small  protrusions  is  several  times  and  up  to  over  ten 
times  the  ideal  thickness  of  the  laminar  flow  boundary  layer.  At 
this  time,  the  roughness  not  only  influences  the  transitions  of 
the  boundary  layer  but  also  influences  the  velocity  and  tempera¬ 
ture  section  structure  of  the  entire  boundary  region.  Further, 
this  can  also  cause  the  heat  flow  on  the  surface  to  increase 
greatly  and  there  can  be  a  several  fold  disparity  with  the  classic 
computation  results. 

Since  the  1970 's,  although  much  experimental  research  has 
been  done  abroad  [2,3,4],  yet  because  the  problems  are  so  complex 
or  because  the  experiments  themselves  are  difficult,  up  to  the 
present  the  inherent  influences  and  mechanisms  are  still  not 
totally  clear.  Because  of  this,  we  measured  several  types  of  bead 
rough  wall  model  heat  flows  with  a  small  scale  hypersonic  wind 
tunnel  to  compare  and  also  measure  the  heat  flow  value  of  the 
smooth  wall  model. 

Because  the  increasing  effects  of  roughness  on  heat  flow  are 
closely  related  to  the  given  Reynolds  regions  as  well  as  the 
geometric  shape  of  the  model,  we  selected  typical  sphere  cone 
shapes.  Five  models  were  used  and  the  rough  bead  diameters  on 
the  surfaces  of  the  models  were  0,0. 3, 0.5, 0.7  and  0.9mm,  the 

Mach  number  is  5,  the  wind  tunnel  total  pressure  p  was  10-45 

2  ^ 
kg/cm  ,  and  the  corresponding  Reynolds  numbers  Rep  were 

(0.8-3. 6)  xlO^.  For  convenience  of  analysis,  we  also  gave  the 
laminar  flow  around  smooth  wall  sphere  cones  under  similar  con¬ 
ditions,  and  the  theoretical  computation  values  of  the  turbulent 
heat  flow. 
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II.  Symbols 


c 

d 

D 

K 

P 

q 

R 

Ren 

Re/ 


s 

T 

t 

u 

P 

S 

0 

M 

<P 


material's  specific  heat  (kilocalories/kg°C) 

rough  bead  diameter  (mm) 

model ' s  base  diameter  (mm) 

equivalent  height  of  roughness  (mm) 

2 

pressure  (kg/cm  ) 

2 

heat  flow  (kilocalories/meter  seconds) 
radius  of  nose  (mm) 

Reynolds  numbers  P#  u^  D /ju# 

local  Reynolds  numbers  on  surface  of  model  P^s /jj,^ 

distance  along  surface  beginning  from  stagnation 
point  (mm) 

Temperature  (K) 

time  (seconds) 

velocity  (meters/second) 

2 

air  flow  or  material  density  (kg/meter  ) 

thickness  of  model  wall  or  boundary  layer  (mm) 

included  angle  between  surface’s  normal  line  and 
model's  symmetrical  axis  (degrees) 

2 

granularity  coefficient  (kg* seconds/meter  ) 
meridian  angle  (degrees) 


Lower  Symbols 

e  outer  boundary  value  of  boundary  layer 

i  different  measuring  point  values 

sq  stagnation  point  value  of  smooth  wall 

t  flow  stagnation  conditions 

00  flow  conditions 


III.  Test  Conditions 

1 .  Wind  Tunnel 

This  experiment  was  carried  out  in  a  second  order  firing 
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free  efflux  type  conventional  hypersonic  wind  tunnel.  In  order 
to  have  a  sufficient  Reynolds  number  variation  range,  we  sel¬ 
ected  relatively  low  experiment  Mach  numbers.  The  concrete 
experiment  conditions  for  Mach  number  5  are  as  follows: 


(!)*  **/■** 

10  1 

20 

30 

45 

R»l>x  10'* 

0.1  I 

l.t 

2.4 

3.6 

T,  K 

473 

Key:  1 


Pt  kg/cm 


Following  the  increases  of  the  total  pressure,  the  Mach 
number  decreased  slightly.  This  can  be  disregarded  in  the 
experiment . 


2.  Experiment  Model 


The  nose  radius  R  of  the  sphere  cone  selected  for  the  exper 
iment  is  27.4  mm,  the  base  diameter  D  is  60  mm  and  to  make  it 
convenient  to  connect  with  the  tail  strut,  the  model's  base  was 
extended  to  the  cylinder  section  (for  details  see  fig.  1) . 


The  models  used  common  heat  transfer  to  test  the  thin  wall 


structure.  The  models  are  electroformed:  on  the  preprocessed 
plastic  (polystyrene)  sphere  cone  model  and  on  uniform  glue 
the  small  plastic  spheres  go  through  a  sieve  with  different 
diameters.  After  electroplating  forms  a  red  copper  matrix  there 
is  nickel-plating  to  corrode  the  copper  film.  This  base  machin¬ 
ing  attains  a  rough  heat  transfer  model  with  a  surface  which  has 
small  spheres  with  different  diameters.  The  small  spheres  are 
stuck  closely  and  firmly  together  and  do  not  allow  piling  up.  The 
glue  cannot  be  too  thick  so  as  to  avoid  having  gaps  blocking  up 
the  small  spheres.  Therefore,  formation  is  relatively  difficult. 
There  are  five  types  of  formed  models;  the  d=0  smooth  wall  model 
and  the  d=0 . 3 , 0 . 5 , 0 . 7  and  0.9mm  rough  wall  models.  The  d  values 
given  here  are  approximate  values.  The  thermocouple  is  welded  on 
the  inner  surface  and  a  spot  welder  is  used  on  the  two  meridian 
surfaces  of  4>  =0°  and  180°.  Altogether  ten  0.2  mm  diameter 
nickel  chromium-nickel  aluminum  thermocouples  are  welded  on.  The 
nominal  thicknesses  of  the  experiment  models  is  1  mm  and  they  are 
limited  by  the  time  of  electroplating.  In  reality,  this  cannot  be 
done,  a  point  which  will  be  further  discussed  below. 

3.  Experimental  Process  and  Instruments 

In  order  to  obtain  very  large  aerodynamic  heating  and  rela¬ 
tively  uniform  initial  temperature  for  the  model ,  during  the 
blowing  process  the  model  is  first  hidden  in  the  stagnant  water 
area  outside  the  "excess  tube"  of  the  wind  tunnel.  Air  is  used 
for  cooling  and  after  the  wind  tunnel  starts  we  close  the  cool 
air.  At  the  time  of  the  experiment,  the  model  is  released  as  a 
free  falling  body  and  because  there  is  a  relatively  large  dispar¬ 
ity  in  inner  and  outer  pressure  during  the  experiment  the  falling 
velocity  is  very  fast.  The  entire  release  throw  is  about  0.2 
seconds  and  the  time  of  passing  from  the  jet  flow  boundary  to  the 
core  flow  is  even  shorter. 
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The  temperature  changes  on  the  surface  of  the  model  are 
sent  from  the  wind  tunnel ' s  examination  equipment  to  the 
Chinese-built  121  computer  for  print  out.  The  examination 
sampling  speed  is  200  points/second,  precision  is  t  0.5°  and 
the  final  data  processing  is  completed  on  the  Chinese-built  320 
computer.  Flow  field  photography  uses  laser  technology  so  that 
one  by  one  conventional  instruments  are  being  phased  out. 

IV.  Data  Processing 

1.  Related  Constants 

After  eliminating  radial  heat  conduction,  the  heat  flow  com¬ 
putation  formula  for  the  common  thin  wall  heat  transfer  model  is: 

Qi-Pcbr^ff-  ( 1 ) 

In  the  formula,  P  and  c  are  the  related  physical  character¬ 
istic  parameters  of  the  model's  material  and  c  is  the  function  of 
the  wall  temperature.  An  approximate  constant  quantity  can  be 
taken  from  our  useful  temperature  range.  Here,  we  take  the  value 
when  the  wall  temperature  is  20°C,  (°  =8900  kg/meter^,  c=0.106 
kilocalories/kg°C . 

b,  is  the  model's  wall  thickness  in  the  thermocoupler's 
welding  point  area.  Because  the  model  has  a  curved  surface  and 
small  spheres,  even  though  the  electrode  shape  makes  a  corres¬ 
ponding  change,  the  wall  thickness  of  the  model  is  still  not 
uniform.  Even  the  smooth  wall  model  also  usually  has  thin  stagna¬ 
tion  points  and  a  thicker  cone  surface.  Furthermore,  the  thermo¬ 
coupler's  welding  point  is  welded  on  the  inner  surface  of  model 
according  to  the  printed  marks  and  the  welding  point.  After  the 
welding  points  are  dropped  on  the  outer  surface  small  spheres  they 
also  fall  in  the  spaces  of  the  small  spheres  which  are  difficult 
to  know.  Therefore,  the  f>,  which  depend  on  direct  measurement 
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or  computation  all  pose  difficulty.  If  the  welding  is  bad,  small 
concave  hollows  may  appear  and  thus  change  the  wall  thickness. 
For  this  reason,  this  experiment  was  calibrated  on  self-made 
static  calibration  equipment  with  a  30  kilowatt  infrared  iodine 
tungsten  lamp  set.  The  standard  for  static  calibration  is  to 
take  a  material  with  known  thickness  of  similar  nickel  plate 
and  a  heat  insulation  plate  is  put  on  the  model  and  nickel  plate 
surface  between  the  commonly  prepared  Chinese  ink,  test  piece 
and  lamp  set.  At  the  time  of  calibration,  the  heat  insulation 
plate  is  quickly  taken  out  and  a  six  pen  recorder  is  used  to 
record  the  temperature' change  curve  of  the  two..  We  can  convert 
the  hi.  of  each  welding  point  on  the  model  from  equation  (1)  . 
See  table  1  for  detailed  numerical  values. 
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Table  1  Coordinates  of  the  Measuring  Points  and  the  Wall  Thick¬ 
ness  of  the  Model  at  Those  Points  (Static  Calibration) 

Key:  1  Sequenae  number  of  measuring  points 
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Key :  2 .  Or 

3.  Sequence  number 

4.  In  the  table,  $  ^  is  the  mean  value  of 

multiple  static  calibration?  the  unit 
is  millimeters 

2.  Computation  of  Derivative  dT^/dt 

Most  often  prior  to  the  experiment  the  temperature  changes 
in  accordance  with  the  time  and  the  shape  of  the  T-t  curve  is 
unknown.  Moreover,  the  curves  differ  in  accordance  with  the 
differences  of  the  experimental  conditions,  model  shape  and  posi¬ 
tions  of  the  measuring  points  of  each  wind  blow.  For  example,  in 
this  experiment,  the  temperature  in  the  vicinity  of  the  smooth 
wall  sphere  cone's  stagnation  points  rises  very  quickly  and  for 
the  most  part  there  are  conic  and  cubic  curve  changes  of  time  t. 
Yet  on  the  cone's  surface,  the  temperature  gradually  rises 
linearly.  It  is  very  difficult  to  find  a  general  equation  nor 
is  it  rational. 

In  consideration  of  the  above  special  points,  we  selected  a 
multiple  regression  analysis  method  commonly  used  in  mathematical 
statistics.  In  principle,  it  can  automatically  select  appropriate 
curve  equations  for  each  measuring  point  from  the  elementary 
function.  Moreover,  the  mean  square  deviation  of  these  curves  is 
minimal.  This  is  different  from  the  method  found  in  reference 
[3]  and  naturally  the  method  is  overelaborate.  This  is  not 
difficult  to  resolve  due  to  checking  and  the  application  of 
computers . 

During  the  experimental  process,  before  the  model  is  in¬ 
serted  into  the  air  flow  the  wall  temperature  is  uniform  and 
after  the  wind  tunnel  starts  there  are  not  very  large  fluctua¬ 
tions.  After  being  inserted,  because  the  sequence  of  each 
measuring  point  of  the  model  making  contact  with  the  air  flow 
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differ  the  times  the  wall  temperature  departs  from  the  initial 
line  are  different.  Furthermore,  some  of  the  heat  flows  of  each 
measuring  points  are  large  and  some  small  and  some  seem  to  be 
invisible.  These  pose  certain  difficulties  for  the  automatic 
processing  of  the  data  and  so  require  proper  formulations.  To 
analyze  the  actual  circumstances  of  this  experiment  and  make  an 
integrated  comparison  we  selected  each  wind  blow  to  be  more  or 
less  close  to  (or  an  average  of)  the  initial  temperature  of  25°C, 
usually  a  little  higher  than  room  temperature.  Afterwards,  the 
last  measuring  point  regression  curve  sticks  in  or  out  to  the 
corresponding  found  time  t  of  the  wall's  temperature  and  is 
fixed  as  the  effective  time  of  the  wall  temperature  beginning 
to  depart  from  the  initial  line.  The  corresponding  dT./dt  is  the 
desired  derivative  value.  The  data  points  in  the  unstable  process 
where  the  T-t  curve's  initial  section  heat  flow  is  too  low  are 
set  aside  beforehand.  Naturally,  they  can  also  be  rejected  in  a 
certain  order  or  corrected  afterwards. 

There  is  a  strong  relationship  between  the  results  of  re¬ 
gression  analysis  and  the  selection  of  the  elementary  function. 

It  is  often  necessary  to  seek  help  from  professional  knowledge 
ana  here  we  selected  a  quintic  multi-item  formula  of  time  t. 

The  T-t  curve  regression  analysis  mean  square  deviation  of  the 
experimental  point  is  about  -  0.5°.  Yet  this  is  only  preliminary 
and  still  awaits  further  improvement.  See  reference  [5,6]  for 
the  mathematical  writing  of  this  method. 

V.  Results  of  Experiment  and  Discussion 

The  results  of  the  experiment,  that  is  the  different  rough¬ 
nesses,  are  given  separately  in  figs.  2-6.  Due  to  the  limit  of 
space,  the  data  table  cannot  be  given  again  and  we  will  only 
quote  from  it  when  necessary.  The  data  in  the  figures  are  the 
regression  computation  results  of  10  experimental  points  obtained 


in  one  second.  For  convenience  of  analysis  and  comparison, 
besides  independently  giving  the  smooth  walls  stagnation  point 
heat  flows,  for  the  rest  we  gave  the  specific  heat  flow 
(q^/qsQ)  forms  equivalent  to  the  smooth  wall's  stagnation  point 
experimental  values. 

1.  Heat  Flow  of  Smooth  Wall 

From  the  form  given  by  the  data,  the  selection  of  qsQ  itself 

is  worthy  of  attention.  Here  we  use  the  mean  value  of  the  cubic 
repeatability  experiment  which  is  relatively  close  to  the 
theoretical  computation  results  for  the  laminar  flow  stationary 
points  given  by  J.A.  Fay  and  F.R.  Riddell  [7] : 


*<**/■*>  fi  j 
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33.45 
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11.10 

27.00 

|  33.10 

j  40.63 

Key:  1.  p  (kg/cm) 

2 

2.  (Kilocalories/meter  seconds) 

3.  Experimental  value 

4.  Theoretical  value 

For  convenience  of  analysis,  we  used  L.  Lees  '  method  [7 ] 
and  the  reference  enthalpy  method  [7,9]  to  make  corresponding 
computations  of  the  partial  laminar  flow  of  the  sphere  cone 
models  and  the  turbulent  flow  and  heat  flow  distribution.  The 
computation  results  are  separately  indicated  with  broken  lines 
in  figures  2-6. 
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Fig.  2  Heat  Flux  Distribution  Along  Surface  for  Smooth  Wall 
Model  With  d=0 


Key:  1.  Theoretical  computation  values 

2.  Turbulent  flow 

3 .  Laminar  flow 

4.  (Kg/cm2) 


Fig.  2  is  a  comparison  of  the  experiment  and  theoretical 
conditions  of  the  distribution  of  the  smooth  wall’s  (d=0)  heat 
flow.  In  the  figure,  the  experimental  points  indicate  the 
results  of  each  experiment  (usually  we  took  the  results  of  the 
first  experiment) .  It  can  be  seen  from  the  figure  that  the  dis¬ 
tribution  of  a  part  of  the  nose's  heat  flow  is  similar  to  the 
computation  values  of  the  laminar  flow.  It  rises  with  the 
enlarging  of  the  wind  tunnel's  total  pressure  p^,  that  is  the 
enlarging  of  the  RPO ,  yet,  in  general,  it  is  still  in  a  laminar 
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flow  state.  This  is  also  very  close  to  the  analysis  of  the  model 

surface’s  partial  Reynolds  number  Re,  computation  results.  When 

p  =45  kg./cm^,  in  the  0=50°,  70°  area  R. u  is  separately 
t  5  5 

5.8x10  and  6.37x10  .  It  can  be  said  that  in  most  situations  it 
is  in  a  transition  boundary  state.  There  is  possible  transition 
on  the  sphere  cone  surface  and  when  in  the  cone  surface's  two 
measuring  points  (measuring  points  numbers  9  and  10)  of  the  T-t 
curve,  a  double  slope  situation  appears.  Theoretical  computations 
could  not  be  made  because  there  were  too  few  measuring  points. 

2.  The  Influence  of  a  Rough  Wall  on  the  Heat  Plow 

To  analyze  the  influence  of  roughness  by  simply  referring 
to  the  heights  of  the  rough  beads  is  not  sufficient.  It  is  also 
related  to  the  distribution  densities  and  configurations  of  the 
beads.  Usually,  the  concept  of  equivalent  roughness  is  used  but 
the  influence  mechanism  is  not  clear  and  it  is  also  difficult  to 
have  a  unified  method.  Here  we  used  the  experimental  results  zf 
reference  [1]  on  the  equivalent  bead  roughness  with  the  aim  that 
it  approximately  comprehends  the  equivalent  measure.  We  also 
take  the  ideal  smooth  wall  laminar  flow  boundary  layer  thickness 
S  as  the  reference  quantity.  When  pt=10  kg/cm  ,  the  computa¬ 
tion  results  are  as  follows: 
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Heat  Flux  Distribution  Along  Surface  For  Rough  Wall 
Model  With  d=0.5  mm 

Key:  1.  Turbulent  flow 

2.  Laminar  flow 
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Fig.  5  Heat  Flux  Distribution  Along  Surface  For  Rough  Wall 
Model  With  d=0.7  mm 

Key:  1.  Turbulent  flow 

2.  Laminar  flow 

3.  (kg/cm^) 
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Fig.  6 


Heat  Flux  Distribution  Along  Surface  For  Rough  Wall 
Model  With  d=0 . 9  mm 


Key: 


1 .  Turbulent  flow 

2.  Laminar  flow 


3.  (kg/cm  ) 


In  figures  3-6,  in  order  for  experimental  points  with  the 
same  roughness  and  different  total  pressure  to  be  put  on  a 
figure,  we  used  logarithmic  coordinates.  Thus,  when  looking 
at  the  figures  one  must  keep  this  in  mind. 


It  is  not  difficult  to  see  from  the  figures  that  under  all  of 
the  experimental  conditions  the  obvious  transition  and  develop¬ 
ment  was  the  turbulent  flow.  After  9  =20°,  all  basically 
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surpassed  their  respective  turbulent  flow  theoretical  values 
(as  shown  in  broken  lines)  and  there  was  a  very  large  increase 
of  the  cone  surface's  heat  flow  as  compared  to  the  results  for 
the  smooth  wall. 

Under  the  same  roughness  conditions,  the  experimental  points 

under  low  total  pressure  were  very  close  to  the  turbulent  flows 

theoretical  broken  line.  For  example,  when  d=0.3,0.5  and  0.7  mm, 

2 

the  results  of  pt=10  kg/cm  were  all  like  this  and  when  d-0.9  mm, 
aside  from  the  stagnation  points,  they  were  also  basically  this 
same  way.  Following  the  enlarging  of  P  and  Re*  the  effects  of 
the  enlargement  of  roughness  on  the  heat  flow  was  even  more 

2 

obvious.  If  d=0.3  mm,  when  total  pressure  pt=10, 20, 30, 45  kg/cm  , 

the  peak  values  of  q./b  were  2.09,2.41,3.09  and  3.44.  The  range 

1  so 

of  increase  was  quite  large. 

Under  different  roughness  conditions  and  similar  total 
pressure,  following  the  enlarging  of  the  roughness  of  the  small 
sphere  diameter  d,  the  heat  flow  value  had  only  a  slight  increase 
or  small  variation.  This  shows  that  the  roughness  of  the  small 
sphere  diameter  is  sufficiently  large.  At  the  same  time,  it  also 
shows  that  even  under  relatively  large  roughness,  because  the 
total  pressure  is  low,  the  roughness  only  affects  boundary 
layer  transition.  Yet,  under  high  total  pressure,  this  type  of 
heat  flow  variation  is  not  large  and  it  possibly  signifies  that 
the  influence  of  roughness  on  heat  flow  increases  carries  certain 
limitations . 

The  rough  wall  heat  flow  distribution  along  the  surface  is 

generally  lower  at  the  two  ends  and  the  peak  value  heat  flow  is 

in  the  vicinity  of  the  sonic  point  region.  Taking  the  d=0.9  mm 

model  as  an  example,  aside  from  the  stagnation  point  the  peak 

value  heat  flow  qVq  in  the  vicinity  of  the  sonic  point  region 

i  so 

approaches  4,  and  when  the  conditions  of  the  stagnation  point 
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differ  from  other  models  including  low  pressure  data,  the  heat 

flow  clearly  increases.  When  the  maximum  total  pressure  is 
2 

p  =45  kg/cm  ,  q./q  approach  6  and  it  seems  that  this  is 

t  1  SO 

possibly  related  to  the  stagnation  point's  local  configuration 
changes . 

Furthermore,  it  can  also  be  seen  from  figures  3-5  that  there 
appears  to  exist  a  transition  region  with  less  than  ideal  turbul¬ 
ent  flow  heating  between  the  stagnation  point  and  0  =20° 
measuring  point.  Whether  or  not  this  signifies  the  existence  of 
a  laminar  flow  region,  the  influence  of  the  roughness  of  this 
region  on  heating  is  not  critical. 

The  influence  of  roughness  on  the  heat  flow  being  this  large 
is  worthy  of  attention  and  naturally  the  problems  are  quite  com¬ 
plex.  Below  we  will  summarize  and  simply  discuss  the  existing 
test  results. 

we  can  see  from  the  data  fluctuations  that  pfc  is  a  parameter 
worthy  of  attention  and  seems  to  cause  an  "enlarging"  effect.  The 
changes  of  pfc  are  also  affected  by  surface  static  pressure  pe 
and  local  Reynolds  number  . 

Under  low  pressure  conditions,  the  roughness  acts  as  the 

laminar  flow  boundary  layer's  unstable  factor  and  this  type  of 

unstable  small  disturbance  strengthens  with  th«~  increase  of 

roughness.  After  the  roughness  is  sufficiently  large,  transition 

appears  under  relatively  low  Res  conditions  and  a  turbulent  flow 

boundary  layer  develops  and  forms  just  as  is  seen  in  the 
2 

Pt=10  kg/cm  data  situation.  The  artificial  transition  rough 
area  method  commonly  used  in  wind  tunnel  experiments  is  based  on 
analagous  principles. 

Following  the  enlarging  of  wind  tunnel  total  pressure  pfc 
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there  is  also  a  corresponding  enlargement  of  model  wall  surface 
static  pressure  pg.  The  small  disturbance  caused  by  the  rough¬ 
ness  seems  to  strengthen  and  this  causes  the  temperature  distribu 
tion  inside  the  entire  surface  boundary  region  to  be  even  more 
uniform.  This  in  turn  causes  an  enlargement  of  the  temperature 
gradient  in  the  vicinity  of  the  surface  which  can  result  in  a 
manifold  increase  of  the  rough  wall  heat  flow.  It  appears  that 
a  more  rational  equivalent  roughness  concept  is  that  the  surface' 
static  pressure  pg  (or  pfc)'  should  be  a  variable  worthy  of  atten¬ 
tion  . 


The  influences  of  Re/  are  easily  understood.  Following  the 
enlargement  of  pfc,  Re*  has  a  corresponding  enlargement  and  the 
thickness  of  the  boundary  layer  then  decreases.  Here,  we  will 
again  use  the  results  of  the  laminar  flow  of  a  smooth  wall.  When 
p  =45  kg/cm2,  &•-••**  0-033  ^  and  thus  the  corresponding  K/r 

has  a  manifold  increase  and  in  essence  enlarges  the  corresponding 
roughness.  Summarizing  these  two  aspects  it  is  not  difficult  to 
conceive  that  the  region  of  the  critical  influence  of  roughness 
on  heat  flow  is  on  the  end  of  the  nose. 
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2 

Fig.  7  Laser  Schlieren  Photoqraph  d=0.9  mm,  pt=20kg/cm 

Looking  at  the  photograph  in  fig.  7  of  the  related  flow  field 
of  the  rough  wall  model  can  help  us  further  discuss  the  problem. 

It  can  be  seen  from  the  photograph  that  there  is  an  obvious 
oblique  shock  waves  group  in  the  vicinity  of  the  sonic  point 
region.  (  B  =30°-40°) .  Their  strength  increases  with  the  enlarge¬ 
ment  of  the  total  pressure  and  roughness.  After  there  is  suffic¬ 
iently  large  roughness  and  waves,  the  rise  of  static  pressure 
Pe  often  produces  complex  shock  wave  boundary  layer  disturbance 
and  separated  reattached  flow  because  of  this.  It  seems  that  this 
phenomenon  is  a  very  important  factor  contributing  towards  the 
many  fold  increase  of  this  region's  heat  flow.  In  the  subsonic 
area  in  front  of  the* sonic  point,  the  roughness  influence  shown 
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in  the  experiment  was  not  an  important  fact  but  seemed  to  be 
very  good  collateral  evidence  for  the  crux  of  the  problem.  The 
transition  under  low  pressure  conditions  mentioned  above  also 
appears  to  be  related  to  the  separated  vorteoc  group  transition 
which  again  proves  the  influence  of  roughness  on  heating  and  the 
configuration.  Naturally,  this  discussion  is  very  preliminary 
and  awaits  further  detailed  experiments. 
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A  GENERALIZED  J- INTEGRAL  OF  TWO  DIMENSIONS 
by  Zhang  Xing  and  Wu  Guoxun 

(Beijing  Institute  of  Aeronautics  and  Astronautics) 

Abstract 

Beginning  from  the  energy  release  rate,  this  paper  derived 
the  generalized  J-integral  for  support  of  surface  forces  and 
body  forces  of  plates  with  variable  thickness  in  non-uniform 
temperature  fields.  It  also  proved  the  conservation  properties 
of  the  generalized  J-integral. 

I.  Preface 

We  know  that  J.R.  Rice  proposed  the  J-integral  for  plates 
with  equal  thickness  under  non-body  force  and  equal  temperature 
conditions  [  1]  and  that  the  J-integral  was  later  extended  and 
applied  widely  in  engineering.  Very  recently,  there  have  been 
articles  which  have  considered  the  influences  of  body  force  and 
temperature  [2,4] . 

Yet,  like  the  importance  of  a  turbine  disk  in  an  aircraft 
engine,  it  not  only  supports  surface  forces  and  body  forces  but 
the  plate  itself  is  a  plate  with  variable  thicknesses  in  a  non- 
uniform  temperature  field. 

Beginning  from  the  energy  release  rate,  this  paper  derives 
the  generalized  J-integral  for  the  support  of  surface  forces 
and  body  forces  of  plates  with  variable  thicknesses  in  non-uniform 
temperature  fields.  It  also  proved  the  conservation  property  of 
the  generalized  J-integral. 
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Due  to  the  limit  of  space,  the  relationship  of  generalized 
J-integral  and  crack  tip  singularity . will  be  discussed  in 
another  article. 

i 

I 

We  have  seen  the  results  of  reference  [2]  from  other  writers 
[3]  but  did  not  see  the  original  article. 

II.  Energy  Expressions  of  Physical  Equations 

1.  Linear  Elasticity 

Based  on  elastic  mechanics,  we  considered  that  the  temperature 
effect  of  the  elastic  stress  and  strain  relation  possesses  the 
following  form: 


+  hi  *  ki  —  tiikiaT 

In  the  above  formula,  CTj is  the  stress  component,  £ _  is  the 

strain  component,  A  and  jU  are  Lame  constants,  k  is  the  volume's 
elastic  constant,  T  is  the  temperature,  a  is  the  linear  expan¬ 
sion  coefficient,  and  is  the  Kronecker  Delta. 


The  unit  volume  of  deformation  potential  energy  W  is  equiva¬ 
lent  to  [5] 


W-§*  ***!'"+- -V  (  2  )  , 

When  formula  (1)  is  substituted  into  formula  (2),  after  integra¬ 
tion  we  can  obtain 

k  3 aTineiH — (  3  ) 
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Based  on  formulas  (1)  and  (3) ,  we  can  prove  that: 


di„ 


At  the  same  time 


alV 

-jf~m »  —  Zkat>„t„+  9kalT 


2.  Elasto-Plastic  Conditions 


(4) 


(5  ) 


Taking  6  and  £  to  separately  represent  the  mean  stress  and 
the  mean  strain,  the  stress  spherical  quantity  and  stress  partial 
quantity  are 

—  /  —  &,/<»  (6) 


and  the  strain  spherical  quantity  and  strain  partial  quantity  are 


tni* 


(7> 


Based  on  HTT^lOUJIlH 
know  that : 


's  total  quantity  theory  [6],  we  can 


k  («</“ 

-f—sH'  I 


'n 


(8) 


In  the  above  formula,  °»  and  **  are  separately  the  equivalent 
stress  and  equivalent  strain.  They  are  equal  to 


0#“/~i t0"0" 

(«) 

M-/  — 

(10) 
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The  deformation  potential  energy  of  the  unit  volume  is  still 
equal  to 


-  (3or)*-J'^K;+oJ/)(Je;f+</e;/)+-|-(3oD*  (11) 


When  we  substitute  formula  (8)  into  formula  (11)  and  take  into 
consideration  formulas  (9)  and  (10) ,  we  can  obtain: 


bnei;bn^iti~  k  3fl7'6/;t(/  +  — | — (3Q71)1 


(12) 


We  can  know  from  formulas  (8)  and  (12)  that: 


aW 


it 


rs 


% 

(13) 


At  the  same  time 


aW 

TjT  ”  -  Zkctbutn+gka'r 


(14) 


I-II.  The  Generalized  J-Integral  and  Energy  Release  Rate 

We  know  that  the  system's  total  potential  energy  TT  is 
equal  to 


n-J  j* c  SiUiids-  ^B,u,tdA  (15) 

In  the  above  formula,  Si  is  the  surface  force,  B<  is  the  body 

force,  ui  is  the  displacement,  A  is  the  area  of  the  plate  and 
Cg  is  the  boundary  of  the  force  (fig.  1) . 
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Fig.  1  Generalized  J-Xntegral  and  Energy  Release  Rate 
Thus,  energy  release  rate  G  is  equal  to 


G 


J_  dU  __J_J  f  o  dm 
i.  da  f.  ( J  c,1  da 


(16> 


In  the  above  formula,  t  is  the  plate  thickness  in  the  crack  tip 
area . 

On  given  displacement  boundary  C  ,  there  is  no  relation 

between  displacements  u,  and  a.  Thus  o .therefore  at  C 

i  da  u 

fc.Srarw'"*  (I7) 

Taking  the  above  formula  into  consideration,  we  can  rewrite  form¬ 
ula  (16)  as  follows: 
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In  the  above  formula,  C  is  the  entire  boundary  of  the  plate. 


We  introduce  a  flow  coordinate  system  (X^X^)  which  takes 

the  crack  tip  as  the  point  of  origin  and  the  crack’s  extended 
line  as  the  horizontal  coordinate.  It  has  the  following  re¬ 
lationship  (fig.  1)  with  the  original  fixed  coordinate  system 

(Xj^ ,  x2)  : 


(19) 


When  compared  to  a  fixed  point,  yy  remain  unchanged.  When  a  has 
changes,  X^  likewise  changes.  Thus,  the  change  of  any  physical 
quantity  P  is: 

dp—£-<ia+jjr;dXi  (2°) 

At  the  same  time,  we  can  know  from  formula  (19)  that: 

dXt--do  (21) 


Thus 


(2  2T 


dP  o  dP  dP 
da  da  dXx 


It  should  be  pointed  out  that  in  the  above  formula, 


dP__jP_  I  dP  _  dP  I 

da  ~  da  I,,.,,  da  ~  da  !*„*, 

dP  dP  I 

zn  =  uxt 


Based  on  the  law  of  composite  function  partial  derivation. 


dP  _  dP  dxx  dP  dXt  dP  dX,  dP 

dXx  dxx  dXx  dXt  dX,  ~  dxt  dX,  ~  dX,  ( 2 3) 


If  formula  (23)  is  substituted  into  formula  (22) ,  we  have 


dP  _  dP  dP 

da  da  dxt  (24) 


Taking  formula  (24)  into  account,  we  can  know  that: 


1 

o 

Ls^ 

(25) 

UA~LBrSrUA 

(26) 

(27) 
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Based  on  the  law  of  composite  function  partial  derivation,  the 
above  formulas  can  be  rewritten  as  follows: 


aw 

A  at,. 


Because  T  is  not  the  function  of  a,  based  on  formula  (24)  we  can 
know  that: 


dT  ^  aT _ *T  _  *T 

da  m  ao  ~  ax,  ~  ’  aa  ~  ax, 


(29) 


If  formula  (29)  is  substituted  into  formula  (28),  we  have: 


c  MildA.c  iSL^L  UA+  c  jl  ua„  r  f,dA  (s0) 

J  a  da  J  a  at,,  aa  j  A  a,  ax,  J  a  ax, 


Taking  this  a  step  further,  based  on  the  partial  integral  form¬ 
ula  and  surface  integral  —  linear  integral  transformation 
formula,  we  can  know  that: 


f  ^L 

J  A  **, 


am 

A  ax. 


*wJkiA 


(31) 


Taking  into  consideration  formula  (31) ,  formula  (30)  can  be 
rewritten  as  follows: 
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(32) 


Based  on  the  principle  of  virtual  work,  we  have 


f  S,  — f<fr+  f  B,-~  id  A*  f  On  atJ'fdA  (33) 

Jc  M  J  A  da  J  ^  da 


Furthermore,  based  on  physical  equation  (4)  or  (13) ,  we  can  know 
that : 


»(/ 


aW_ 


(34) 


Thus 


r  Sr^-tds+  r  B<~i~tdA-  r  —  (35) 

J  c  da  J  A  da  J  A  to,,  da 

If  formulas  (25),  (26)  and  (32)  are  substituted  into  formula  (18) 

and  we  take  into  consideration  formula  (35) ,  we  have: 
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(36) 


G - J~\  \  Wtdxt-  f  W  p-dA-  f  -j^tdA 

t.  U  C  J  ,4  <»*,  J  ^  <»7  iXy 

-  f  sr?-  id*~  f  Bc~r-idA ) 

JC  «i  J  A  <>*,  ) 


The  right  side  of  the  equal  sign  is  the  generalized 
J-integral  when  the  path  of  integration  is  the  boundary  of  the 
plate. 


IV.  The  Definition  of  the  Generalized  J-Integral  and  its  Path 
Independent  Property 


The  r  curve  is  taken  as  a  curve  which  begins  from  the  lower 

surface  of  the  crack  and  ends  on  the  upper  surface  of  the  crack 
including  the  crack  tip.  Letting  ££  represent  the  surrounding 
area  of  the  F  curve  and  upper  and  lower  surfaces  of  the  crack 
(fig.  2)  then  the  J-integral  is  defined  as  follows: 
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Fig.  2  Generalized  J-Integral  and  Its  Path  Independent 
Property 

Below  we  will  prove  that  there  is  no  relation  between  the 

generalized  J-integral  and  its  path.  For  this  reason,  the  p* 

curve  is  a  curve  which  begins  from  the  lower  surface  of  the 

crack  and  ends  on  the  upper  surface  of  the  crack  including  the 

crack  tip.  Q  represents  the  surrounding  area  of  the  p' 

curve  and  upper  and  lower  surfaces  of  the  crack.  We  let  P 

represent  the  closed  curve  formed  from  T  ,  V  as  well  as  the 

upper  and  lower  surfaces  of  the  crack,  and  Q  *  represent  the 

* 

surrounding  area  of  |  .  When  we  consider  that  on  the  upper 

and  lower  surfaces  of  the  crack 

(38) 

0  (39) 

then  based  on  the  definition  of  the  J-integral,  we  have: 
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(40) 


Based  on  the  linear  integral — surface  integral  transformation 
formula,  we  have: 


Sr***’  JV^T 


(41) 


Based  on  the  law  of  composite  function  partial  derivation. 


Wi) 


aW 

at,, 


■**/; 

*xt 


+ 


•T  +r  " 


ax 


ax. 


C42) 


If  formula  (42)  is  substituted  into  formula  (41) ,  we  can  know 
that : 


j >"*•  -  f  -Sr  r „•  %  £- <ja+ f  ^ 


<«> 


If  formula  (43)  is  substituted  into  formula  (40) ,  we  can  obtain 

Based  on  physical  equation  (4)  or  (13) ,  we  can  know  that: 


*W 


-O/l 


(45) 


Further,  based  on  the  principle  of  virtual  work,  we  have: 
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T - r — » 


(48) 


If  formula  (45)  is  substituted  into  formula  (46)  and  then  form¬ 
ula  (46)  is  substituted  into  formula  (44),  we  can  obtain: 

y-J'-O'  (47) 

If  the  attained  surface  force  and  body  force  for  a  plate 
with  variable  thickness  is  located  in  a  non-uniform  temperature 
field,  the  generalized  J-integral  also  has  no  relation  to  the 
path  of  integration. 


V.  Application  Examples  of  Generalized  J-Integral 


After  applying  the  conventional  finite  element  method  to 
find  the  nodal  point  displacement,  we  can  then  calculate  the 
value  of  the  generalized  J-integral.  We  know  that  under  linear 
elastic  conditions  the  energy  release  rate  is  equal  to 

-  •  ■ 

.  G— frOCl+XI)  (¥»£#) 5  (48) 

Key:  1.  Plane  Stress 


In  the  above  formula, 
Because 

therefore 


and  Kj  are  the  stress  intensity  factors. 


G  =  J 


(49) 


/  — ^(tf!+/n) 


(60) 
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If  the  crack  is  located  on  the  component's  geometrically 
symmetrical  line,  we  can  use  the  principle  of  superposition  to 
find  Kj  and  K^;  when  the  crack  is  located  on  a  non-geometrically 
symmetrical  line,  we  can  find  the  stress  intensity  factors  under 
non-composite  conditions. 

VI.  Conclusions 

1)  The  generalized  J-integral  can  be  defined  as: 


■4-{J 

-J> Sr'^-L 


2)  According  to  the  above  definition,  the  J-integral  has  the 
property  of  path-independence. 

3)  Based  on  the  above  property,  the  J-integral  is  equal  to 
the  energy  release  rate  of  the  system. 

4)  in  consideration  of  the  above  property,  we  can  calculate 
the  stress  intensity  factors  after  the  nodal  point  displacements 
have  been  obtained  by  the  finite  element  method  in  linear 
elastic  fracture  mechanics. 

5)  In  elasto-plastic  fracture  mechanics,  we  obtained  results 
similar  to  those  of  Huchinson  [7],  that  is  the  stress  field  in 
the  vicinity  of  the  crack  tip  was  determined  by  the  generalized 
J-integral.  Thus,  the  generalized  J-integral  can  be  used  to 
establish  a  fracture  criterion  under  elasto-plastic  conditions. 
This  problem  was  prepared  according  to  precedent  and  it  is 
discussed  in  another  article. 


59 


References 


[  1]  H.Liebowitz: Fracture  Vol.  I  P.  192-308  1968. 

[2]  Ban  Tian:  See  references  mentioned  in  reference  [3] . 

[3]  Gong  iBenbo,  Tufan  Minggong,  Jigang  Zhunfu  et.al.  The 
Cracking  Strength  of  Rotating  Spherical  Plates,  Collection 
of  Essays  from  the  Mechanics  Meetings  of  Japan  (A  edition) , 
Vol.  45,  No.  396  (297u-82) . 

[4]  W.K.  Wilson  and  L.W.  Yu:  The  use  of  J-integral  in  thermal 
stress  crack  problems.  International  Journal  of  Fracture 
Vol.  15  No.  4  1979. 

[5]  Boley  Weiner:  Theory  of  thermal  stress  P.262-269  1960. 

[6]  A.  A.  Hitnui^nacnagm  P. 97-N  19«t. 

[7]  J.w.  Huchinson:  Singular  behaviour  at  the  end  of  a  tensil 
crack  in  a  hardening  material.  J.  Mech.  and  Phys.  of  Solids. 
Vol.  16  P.13-18  3968. 


60 


DESIGN  CONSIDERATIONS  ON  GLASS  FIBER  REINFORCED  PLASTIC 
COMPOSITE  STRUCTURES 

by  Zhu  Yiling 

(Shanghai  GRP  Research  Institute) 

Abstract 

when  using  glass  fiber  reinforced  plastic/contposite  materials 
for  structural  design,  it  is  at  least  necessary  to  take  the  fol¬ 
lowing  three  points  into  account:  designability  of  the  material's 
properties;  non-restriction  of  the  structure's  configuration;  and 
the  linear  dispersion  strength  of  the  material. 

Composite  materials  possess  many  excellent  properties,  yet 
they  are  not  "only  of  one  family."  They  must  also  exist  and 
develop  in  competition  with  other  materials.  If  they  can  simul¬ 
taneously  us-3  two  or  more  of  the  above  strong  points,  their 
future  is  bright.  The  commonly  known  glass  fiber  reinforced 
plastic  aircraft  nose  radar  cap  used  the  strong  points  of  compos¬ 
ite  materials  transmitting  radar  waves  and  high  specific  strength, 
and  the  cap  has  not  deteriorated  after  40  years.  There  has  been 
vigorous  research  carried  out  in  the  field  of  composite  mater¬ 
ials,  especially  in  the  areas  of  varieties  of  raw  materials  and 
their  applications.  There  have  also  been  great  developments  in 
researching  properties,  yet  proper  attention  has  not  been  given 
to  the  long  term  accumulation  of  reliable  design  data.  There 
are  many  structures  which  are  made  of  composite  materials  most 
of  which  are  still  in  the  "replacement  design"  stage,  that  is, 
under  conditions  where  the  shape,  load  and  application  environ¬ 
ment  of  the  prototype  materials  is  unchanged,  it  is  necessary 
to  make  full  use  of  the  strong  points  of  composite  materials 


when  using  composite  materials  to  replace  prototype  materials 
when  manufacturing  parts. 

This  paper  only  takes  into  account  that  when  designing  a 
structure  with  composite  materials,  it  is  at  least  necessary  to 
pay  attention  to  the  following  three  points:  designability  of 
the  material's  properties;  non-restriction  of  the  structure's 
shape;  linear  elasticity  and  the  dispersion  strength  of  the 
material . 

I.  Designability  of  Material's  Properties, 

The  designability  if  a  composite  material's  properties  re¬ 
lates  to  the  forming  of  two  large  categories  of  materials  for 
composite  materials,  basic  and  reinforcement  materials,  which 
can  be  selected  according  to  design  requirements.  It  is  necessary 
that  the  properties  of  both  the  basic  and  reinforcement  materials 
be  understood  clearly  beforehand.  The  selection  of  different  re¬ 
inforcement  and  basic  materials  as  well  as  their  content  ratio 
and  various  layer  forms  can  make  up  various  composite  materials 
with  different  properties.  Taking  the  ratio  of  shearing  elastic 
modulus  G  and  longitudinal  elastic  modulus  E  as  an  example,  if 
we  select  conventional  materials,  then  G/E  is  a  constant  for  the 
specially  designated  materials.  If  we  select  and  use  fibrous 
composite  materials  then  G/E  is  a  designable  numerical  value 
whichigives  designers  an  extra  degree  of  freedom.  Naturally,  the 
various  property  values  of  the  composite  materials  are  conditioned 
by  the  selected  materials. 

The  composite  mechanism  of  these  two  large  categories  of 
materials  are  basically  physical  composites  and  thus  the  proper¬ 
ties  of  composite  materials  can  be  estimated  according  to  the 
properties  of  the  component  materials.  This  can  generally  be  sim¬ 
ulated  by  using  the  parallel  type  model  and  series  type  model  as 
shown  in  fig.  1. 
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Fig.  1  Primary  Models 


Key:  1.  Basic  body 

2.  Fiber 

3.  Parallel  type 

4 .  Fiber 

5.  Basic  body 

6.  Series  type 


For  example,  the  elastic  modulus  along  the  fiber  direction  is 
calculated  by  the  parallel  type  model: 

E  =E -V.+E  V, 
c  f  f  m  f 

In  the  formula,  E  represents  the  elastic  modulus,  V  indicates 
the  volume  content,  and  the  lower  symbols  c,f  and  m  successively 
indicate  the  composite  material,  fiber  and  basic  body  (all  of 
which  will  be  similar  hereafter) .  They  are  very  similar  to  the 
test  values  and  the  longitudinal  strength  can  also  be  estimated 
by  analagous  formulas.  The  elastic  modulus  of  the  vertical  fiber 
direction  can  be  estimated  by  the  series  type  model: 


1_ 

E 


m 


m 
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This  estimated  value  is  lower  than  the  measured  data  and  so 
to  more  accurately  reflect  the  real  conditions  we  can  further 
use  the  series  parallel  composite  model  [1].  Other  estimations 
such  as  for  heat  and  electric  properties  can  also  use  the  same 
type  of  model.  This  type  of  composite  takes  the  properties  of  an 
independent  component  material  as  the  criterion  and  does  not  con¬ 
sider  their  boundary  surface  influences.  However,  some  properties 
focus  on  boundary  surface  effects  such  as  fracture  toughness  but 
the  above  mentioned  series  parallel  model  cannot  be  used  for 
simulation . 

In  consideration  of  application  requirements,  composite 
materials  can  be  designed  into  three  typical  forms  [2] : 

1)  Isotropic  Composite  Material.  It  can  be  composed  of  short 
cut  fibers  and  also  use  0°,  -  60°  or  0°,  -  45°  and  90°  direction 

continuous  fibers.  There  are  also  3,4  and  7  direction  composit- 
materials.  The  strength  and  stiffness  of  this  type  of  material 
is  relatively  low  and  it  is  suitable  for  use  on  structures  with 
relatively  large  load  randomness. 

2)  Orthogonal  Balanced  Composite  Material.  It  can  be  composed 
of  0°  and  90°  direction  continuous  fibers  and  can  also  use  inter¬ 
woven  fabric  for  formation.  There  is  relatively  high  strength 
and  stiffness  in  the  two  main  directions  (longitudinal  and 

latitudinal  directions) ,  yet  the  strength  and  stiffness  of  the 
45°  angle  it  forms  with  the  main  direction  is  relatively  weak. 

It  is  suitable  for  structures  with  strength  and  stiffness  re¬ 
quirements  in  the  two  vertical  directions,  for  example  the  plate 
and  shell.  We  can  still  use  this  type  of  composite  material  for 
plates  (shells)  with  different  frame  and  rib  distances  yet  its 
main  direction  is  not  suitable  to  be  parallel  with  the  frame  and 
rib  but  should  form  a  45°  with  it.  In  this  way,  when  in  a  large 
degree  of  deflection  it  is  even  more  advantageous  to  bringing 
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out  the  strength  properties  of  the  fibers.  It  also  has  a  strong 
shearing  strength  in  the  plate  (shell)  surface.  China's  design 
of  a  wooden  bed  frame  strung  with  crisscross  coir  ropes  is  a 
good  example  of  the  use  of  fibers. 

3)  Unidirectional  Composite  Material.  This  is  a  type  of  com¬ 
posite  material  wherein  the  entire  fiber  is  placed  in  the  same 
direction.  There  is  especially  high  strength  and  stiffness  in 
the  fiber  direction  yet  it  is  weakest  in  the  vertical  fiber 
direction.  We  can  only  use  this  type  of  composite  material  when 
the  load  conditions  are  very  clear,  for  example  for  cable.  In 
most  calculations  beams  and  columns  have  simple  direction  stress 
analysis.  This  is  because  the  second  direction  strength  of 
isotropic  materials  is  not  sufficiently  large  and  thus  it  is  not 
necessary  to  calculate  the  stress  of  the  second  direction.  Taking 
a  square  tube  sectional  beam  of  a  wall  with  equal  thickness  as 
an  example,  we  take  unit  beam  length  ds  as  the  detached  body  (as 
shown  in  fig.  2) ,  the  bearing  moment  of  bending  on  the  section  is 
M,  bearing  pressure  stress  Cf  is  on  the  upper  surface  plate, 

and  the  bearing  pressure  stress  is  on  the  lower  surface  plate. 

Because  the  beam's  bend  and  the  bearing  pressure  on  the  upper 
surface  plate  are  not  on  one  straight  line,  the  pulling  force  on 
the  lower  surface  plate  is  also  on  different  lines. 


Fig.  2  Square  Tube  Subjected  to  Bending 
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They  form  the  pressure  of  a  pair  of  opposite  pressurized  tubes. 
Their  degree  of  accumulation  q  is: 


7  -  JgiLf  \  _o,i  XM\  .a(P;/,)/(6/2)  2o\(  t  \ 

7  dz  Up;  p  e,j,  a‘,=  'EJT~  a,t 

The  largest  bend  M2  in  the  square  tube  section  is  obtained  by 
using  any  structural  mechanics  method 


if** 


The  largest  stress  is 


6  Mt 
'  r 


6 

16 


The  calculations  of  the  upper  surface  show  that  the  beam,  column 
and  rod  also  only  require  second  direction  fiber.  When  estimating, 
not  calculating  the  basic  body's  strength  and  hardness,  we  ob¬ 
tained  the  strength  and  hardness  of  the  composite  material  as 


n„ 


El*  n,+n, 


Jo, 


d  + 

18 


nl+n2  d°' 


=  - 

B  n. +n, 

X  4m 


In  the  formulas. 


n^  and  n2  are  the  first  and  second  direction  fiber  quantities  of 
the  composite  material .  E  and  represent  the  elastic  modulus 
and  strength  of  the  fibers.  From 


we  obtained  the  ratio  of  the  two  directions  of  fibers 


«i  0,  4  E7\  t  ) 

Other  sectional  forms  can  also  use  this  method  for  processing. 
Only  if  there  is  a  rationally  disposed  fiber  quantity  can  the 
strength  properties  of  the  fibers  be  brought  into  full  play. 
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The  strength  and  elastic  modulus  of  the  fibers  are  far 
greater  than  the  strength  and  elastic  modulus  of  the  basic  body, 
that  is 

^  mB  and  Ef  »  E^.  When  in  the  preliminary 

design,  by  not  calculating  the  strength  and  elastic  modulus  of 
the  basic  body  excessive,  errors  will  not  occur.  In  this  way, 
the  strength  design  of  the  structure  can  be  summarized  in  one 
sentence:  arrange  a  fixed  quantity  of  fibers  in  a  certain  direc¬ 
tion.  This  is  sometimes  called  network  analysis.  When  linked  to 
technical  feasibility  and  load  randomness,  it  is  not  possible  to 
place  all  of  the  fibers  on  the  main  stress  direction  and  it  is 
only  possible  to  place  the  fibers  in  the  main  sttess  direction 
within  a  certain  area.  The  stress  in  this  area  can  still  have 
Cf  ^ ,  C?  an<^  ^*12'  anc^  besides  requi  ring  fibers  in  one  and  two 
directions  it  is  also  necessary  for  the  45°  direction  fibers  to 
have  bearing  shearing  force.  The  example  shown  in  fig.  3  illus¬ 
trates  this. 


Fig.  3  Stress  Ratio  Vs.  Fiber  Volume  Ratio 
The  Stress  ratio  is  :  C£> : 

®  :  0  •  2 '  the  fiber 

volume  ratio  is  1:0.5:  (0.2+0. 2) 

After  first  selecting  the  fiber  volume  ratio,  the  formation  of 
this  layered  plate  is  known  and  afterwards  the  stress  of  a  single 
layer  is  calculated  according  to  the  analysis  method  for  layered 
plates.  Further,  a  suitable  strength  criterion  is  selected  for 
examinations  and  revisions  are  carried  out  so  that  thickness 
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reaches  the  thinnest  state.  The  properties  of  composite  layer 
plates  are  superior  to  the  properties  of  single  layer  plates. 

Attention  should  be  given  so  that  the  resin  base  composite 
material  has  three  low  property  values:  tensile  strength  (5 

T  B 

of  the  vertical  fiber  direction;  and  shearing  strength 

ltb  and  Shearing  elastic  modulus  G^T  of  the  vertical  fiber 
direction.  The  use  of  unidirectional  thread  composite  layered 
plates  is  not  necessarily  the  more  appropriate  plan.  The  use 
of  fabric  with  the  same  longitudinal  and  latidudinal  directions 
or  (and)  fabric  with  different  longitudinal  and  latitudinal 
directions  as  the  reinforcement  material  is  sometimes  suitable. 

II.  Non-Restriction  of  Structure's  Configuration 

At  present,  the  structural  designs  of  glass  fiber  reinforced 
plastic  composite  materials  are  still  in  the  "replacement  design" 
stage.  During  the  preliminary  stage,  this  type  of  design  is 
the  only  way.  Yet  it  should  be  stated  that  it  is  not  a  good 
method  owing  to  the  fact  that  when  there  is  replacement  of  indiv¬ 
idual  parts,  they  are  constrained  by  the  surrounding  structures 
so  that  full  play  cannot  be  given  to  the  points  of  the  composite 
materials  for  the  simultaneous  design  of  materials  and  structure. 

Prior  to  molding,  the  fibers  are  soft,  the  basic  body  is 
flowing  and  the  supporting  mold  forms  the  composite  material 
product.  It  is  unlike  the  customarily  used  metal  or  wood  mater¬ 
ials  which  require  the  assemblage  of  complete  materials.  Because 
of  this,  the  size  and  thickness  of  the  component  is  not  influenced 
by  becoming  full  sized  but  can  be  determined  according  to  design 
requirements.  When  beginning  to  consider  structural  designs  on 
composite  materials,  we  should,  as  far  as  possible,  join  the 
joinable  structural  components  into  a  part.  For  example,  beam, 
rib  and  plate  structures  do  not  necessarily  need  to  first  mold 
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beams,  ribs  and  plates  but  later  we  can  carry  out  the  putting 
together  so  that  one  composite  part  is  formed  at  a  time.  In 
this  way,  we  decreased  the  joining,  the  amount  of  assemblage 
work  and  the  mold  lightened  the  weight  of  the  structure,  and 
also  raised  the  quality  of  the  product,  thus  lowering  the  cost. 
Fig.  4  is  a  full  GRP  antenna.  It  is  composed  of  6  ribbed  fan 
surfaces,  one  central  circular  surface  and  one  joined  cone.  The 
prototype  structure  is  a  steel  structure  composed  of  a  central 
circular  hole,  radiating  beams  and  circular  rings  which  are  in¬ 
stalled  on  a  GRP  honeycomb  sandwich  construction  reflecting 
surface  or  metal  reflecting  surface.  This  type  of  structure  is 
relatively  light  in  weight  and  the  GRP  only  acts  as  the  base 
material  of  the  metal  network.  It  does  not  bring  the  special 
features  of  the  GRP  lightness  and  high  strength  into  full  play. 
The  modified  design  assimilated  the  structural  layout  scheme  of 
the  metal  structure  and  used  the  GRP's  special  features  of 
lightness,  high  strength  and  easy  formation.  This  decreased  the 
number  of  single  structural  components  and  the  amount  of  assembly 
work.  At  the  same  time,  this  caused  the  weight  of  each  component 
to  decrease  so  that  they  can  be  lifted  by  a  person  and  put  on  the 
back  of  a  horse  .  If  the  strong  points  of  the  rotating  shell  are 
brought  into  full  play,  this  can  cause  the  structure  to  be  even 
more  succinct.  Furthermore,  if  a  helicopter's  rotary  wing  uses 
composite  materials,  it  is  easy  to  manufacture  a  twisted  rotor 
blade  to  replace  the  conventional  level  and  straight  rotor  blade 
and  this  raises  aerodynamic  performance  [3] . 
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Fig.  4  Full  GRP  Antenna 

The  stiffness  of  the  structure  is  not  only  related  to  the 
elastic  modulus  of  the  material  but  also  to  the  configuration  of 
the  structure.  To  raise  stiffness  it  is  necessary  to  consider 
two  aspects.  The  first  was  already  given  attention,  the  finding 
many  types  of  high  modulus  materials?  the  second  has  not  been 
given  enough  attention.  It  was  mentioned  previously  that  glass 
fiber  reinforced  plastic/composite  materials  have  three  low 
material  property  values  and  although  we  can  avoid  using  these 
low  values  in  lamination  designs,  should  rather  use  the  relatively 
high  values  and  yet  it  takes  as  its  cost  the  sacrifice 

of  the  longitudinal  elastic  modulus  and  strength.  When  necessary, 
we  can  also  use  this  type  of  processing  so  that  it  slightly 
loses  some  longitudinal  properties  in  exchange  for  raising  the 
shearing  properties.  By  changing  the  structural  configuration 
we  can  possibly  attain  better  advantage,  be  able  to  manufacture  a 
closed  type  rod  but  not  an  open  type  rod,  and  be  able  to  select 
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a  thick  type  not  a  thin  type  structure.  Taking  a  diving  board 
as  an  example,  a  metal  diving  board  is  a  comb  type  while  one 
made  of  wood  is  a  solid  board.  When  designing  a  glass  fiber  re¬ 
inforced  plastic  diving  board,  if  we  use  a  wooden  board  proto¬ 
type,  then  the  weight  is  excessive?  if  we  use  a  metal  structure 
prototype  then  very  large  twisting  and  swaying  can  occur.  The 
square  tube  section  ( Lx±.i.j.i3  )  is  much  more  suitable.  To  give 
another  example,  because  considerations  for  the  initial  struct¬ 
ural  configurations  of  the  aircraft  nose's  radar  cap  were  not 
complete,  it  was  not  able  to  pass  static  water  pressure  tests. 
Afterwards,  although  it  was  strengthened  it  still  lost  electronic 
properties.  If  considerations  were  more  complete  during  prelim¬ 
inary  designs  each  radius  of  the  cap's  curvature  would  be 
suitably  adjusted  and  at  the  same  time  we  could  completely 
satisfy  the  mechanical  and  electronic  requirements. 

The  properties  of  fiber  composite  materials  possess  very 
strong  directionality  and  even  if  configuration  selection  is 
correct,  if  at  the  time  of  formation  we  do  not  control  the  fiber 
direction  there  can  be  a  loss  of  effectiveness.  There  can  also 
be  difficulty  in  analyzing  strain  data  during  the  structural 
tests  because  of  the  unknown  elastic  modulus  of  the  measuring 
points.  Taking  again  the  above  mentioned  aircraft  nose  cap  as 
an  example,  to  fully  utilize  the  properties  of  the  fibers  we 
separated  five  regions  during  the  formation  of  the  cap  and 
among  four  of  these  five  we  made  every  effort  to  use  the  coin¬ 
cidence  of  the  fiber  direction  and  buckling  axis  to  raise 
buckling  resistance  capabilities.  For  the  other  region  we  made 
the  isotropy  accurate. 

III.  The  Linear  Elasticity  and  Strength  Dispersion  of  the 
Materials. 

The  stress-strain  relation  up  to  fractures  of  many  fibers 
and  other  composite  materials  are  linear  and  only  in  the 
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unfractured  range  is  there  very  weak  "plasticity."  This  is  caused 
by  the  gradual  fracturing  of  the  fibers.  It  is  different  from 
metal  and  therefore  this  type  of  composite  material  has  no  ap¬ 
parent  plasticity  area.  Thus  it  does  not  have  the  capabilities 
to  alleviate  the  defective  area's  high  stress  so  that  relatively 
large  strength  dispersion  occurs.  People  have  frequently  blamed 
composite  material  properties  as  being  unstable  and  when  a  com¬ 
posite  material  breaks  down  quite  a  few  material  properties  are 
unstable  yet  rarely  do  people  consider  the  problems  of  design 
or  material  properties  selection  indices.  Actually,  whether  a 
material  's  properties  are  stable  or  not  is  related  to  the 

selection  index.  Taking  the  yield  limitations  of  number  3  steel 

2 

as  an  example,  the  smallest  value  measured  is  17  kg/m m  ,  maxmimum 

value  is  42  kg/rran^  at  a  difference  of  2^  times.  If  we  take  30  kg/mm 

as  the  yield  index  for  nuirber  3  steel,  it  is  feared  that  many  will 

say  that  the  properties  of  the  steel  are  unstable.  However, 

2 

the  use  of  24  kg/mm  as  the  yield  index  guarantees  an  over  99 
percent  probability  of  exceeding  this  index  and  thus  people  will 
say  the  properties  of  the  steel  are  stable.  This  type  of  index 
should  also  be  fixed  for  composite  materials.  Taking  non  (low)- 
alkali  glass  fiber  as  an  example,  when  its  strength  discrete 
coefficient  is  about  10  percent,  if  we  strictly  control  the  man¬ 
ufacturing  techniques  of  the  composite  materials,  its  discrete 
coefficient  can  be  smaller  than  10  percent.  This  is  because  the 
probability  of  the  weakest  point  of  each  fiber  appearing  simultan¬ 
eously  on  a  section  is  very  small.  However,  because  in  most  glass 
fiber  reinforced  plastics  the  partial  difference  of  the  fiber 
direction  and  the  resin  content  are  not  uniform,  this  causes  the 
discrete  coefficient  to  reach  15  percent  and  when  compared  to 
the  8  percent  discrete  coefficient  of  the  number  3  steel's  yield 
limitation,  the  property  dispension  of  the  glass  fiber  reinforced 
plastic  materials  is  a  little  greater.  It  is  necessary  for 
structural  designers  to  recognize  that  this  dispersion  is  an 
objectively  existing  and  suitably  correct  limiting  performance 
value.  The  two  limiting  values  for  criterion  in  the  standard  5A 


used  by  the  United  States  military  are:  one  is  the  "A"  criterion 
which  takes  a  performance  value  guaranteed  99  percent  probability 
under  95  percent  believability  as  the  limiting  value.  Generally 
used  after  efficiency  is  lost  it  causes  a  situation  wherein  the 
totality  of  the  structure  is  lost.  The  other  is  the  "B"  criter¬ 
ion  which  takes  a  performance  value  guaranteed  to  exceed  90  per¬ 
cent  probability  under  95  percent  believability.  Generally  used 
after  efficiency  is  lost  the  outer  load  is  redistributed  to  other 
bearing  force  components  [4] .  These  types  of  considerations  are 
also  suitable  for  composite  materials.  After  determining  the  lim¬ 
iting  value,  it  is  also  necessary  to  formulate  a  safety  factor. 
British  standards  consider  six  factors.  The  first  is  basic 
safety  factor  3  and  the  other  five  are  factor  1.4-3  related  to 
production  methods,  1.2-2  related  to  bearing  time,  1-1.25  related 
to  operating  temperature,  1.1-2  related  to  alternating  frequency, 
and  factor  1. 1-1.5  related  to  the  solidification  system.  The 
product  results  range  from  6  to  40  [5] .  The  use  of  these  types 
of  safety  factors  cause  the  strong  points  of  the  composite  mat¬ 
erials,  lightness  and  high  strength  not  to  be  manifested..  In 
the  NACA  Lewis  Research  Center  of  the  Structural  Composite 
Materials  Industrial  Company  (SCI  Inc.)  of  the  United  States,  the 

operating  pressure  of  the  modelled  epoxy  Kevlar  49  spherical 

2  2 
container  is  1.9  kg/mm  ,  the  blast  pressure  is  3.27  kg/mm  [6], 

and  the  corresponding  safety  factor  is  1.7.  When  the  glass  fiber 
reinforced  plastic  diving  board  was  mentioned  previously  for  the 
preliminary  design,  we  first  used  safety  factor  10  for  the  dead 
load  design  and  after  the  dynamic  load  was  checked  and  we  ob¬ 
tained  a  safety  factor  of  2.37,  the  selected  safety  factor 
appeared  a  little  smaller  when  fractures  appeared  in  the  small 
batch  of  plates  used.  It  should  be  determined  which  composite 

material  safety  factors  are  worth  examining.  Taking  glass  fiber 
reinforced  fiber  as  an  example,  even  though  we  have  a  fixed 
quantity  of  products,  nevertheless  most  of  the  designs  are 
reasonable.  Some  don't  even  have  the  properties  of  raw  materials,, 
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and  the  lack  of  daily  records  causes  there  to  be  no  way  of 
carrying  out  analysis  after  structural  efficiency  is  lost. 
Designers  should  know  in  relative  detail  the  load  and  utiliza¬ 
tion  environment  and  should  also  obtain  relatively  reliable 
material  property  data  from  sufficiently  large  subsamples  so  as 
to  formulate  relatively  reasonable  safety  factors. 

We  already  know  that  fiber  strength  is  related  to  the  defects 
which  exist  in  the  fibers  and  the  distribution  of  the  defects 
are  random.  The  longer  the  fiber  length,  the  greater  the  prob¬ 
ability  of  defects  appearing  and  thus  its' strength  will  be  lower. 
On  the  contrary,  when  fiber  length  is  shorter  there  is  a  smaller 
probability  of  large  defects  appearing  and  the  strength  can  then 
be  higher.  It  is  very  natural  that  when  designing  composite  mat¬ 
erial  structures  we  must  consider  the  relationship  of  the 
structure's  measurements  and  the  material's  strength.  Taking  the 
example  shown  in  fig.  5,  three  simply  supported  beams  are  sub¬ 
jected  to  maximum  bending  moment  j  PL  and  we  can  no  doubt  select 
equal  measurement  beams  among  conventional  materials. 


Fig.  5  Three  Simply  Supported  Beams  Subject  to  Equal  Bending 
Moment 

This  can  also  be  done  for  composite  materials.  The  most  suitable 
method  should  consider  thesize  of  the  area  occupied  by  an  equal 
stress  level  and  the  largest  stress  area  of  the  three  points  of 
the  additional  load  beam  are  much  smaller  than  the  largest  area 
of  the  four  points  of  the  additional  load  beam.  Secondly,  the 
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large  stress  area  is  also  smaller  than  the  evenly  distributed 
load  beam.  Thus,  the  sectional  measurements  of  these  three 
beams  should  be  different  and  the  measurements  of  the  latter  two 
should  be  larger  than  the  three  points  of  the  additional  load 
beam.  These  problems  are  not  very  outstanding  on  small  measure¬ 
ment  structures,  yet  consideration  should  be  given  on  large 
scale  structures.  When  handling  these  problems,  we  can  use  the 
weakest  chain  analytical  method  to  establish  the  relationship  of 
strength  and  length  as  shown  in  fig.  6. 


Fig..  6  The  Relationship  of  Strength  and  Length 

Key:  1.  Standard  scale  distance 
2 .  Standard  strength 

This  relationship  is  also  suitable  for  width.  For  the  relation¬ 
ship  of  strength  and  thickness  we  can  consider  that  the  proba¬ 
bility  of  the  weakest  position  of  each  thin  layer  concentrating 
on  the  same  section  is  very  small.  When  thickness  increases,  the 
strength  can  be  raised  close  to  the  mean  strength  as  shown  in 
fig.  7  and  the  use  of  the  weakest  strength  of  the  thin  layer  for 
design  tends  to  be  safe. 
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Fig.  7  The  Relationship  of  Strength  and  Thickness 
IV.  Conclusion 


To  design  well  a  composite  material  structure,  it  is  first 
necessary  to  understand  the  outer  loads  clearly,  especially  the 
secondary  loads.  This  is  because  the  second  direction  strength 
of  single  direction  composite  materials  is  very  weak.  Further¬ 
more,  it  is  also  necessary  to  sufficiently  understand  the  oper¬ 
ational  environment  of  the  structure  so  as  to  be  able  to  select 
suitable  base  body  and  reinforcement  materials.  Secondly,  it  is 
also  necessary  to  have  full  knowledge  of  the  material's  proper¬ 
ties  as  well  as  needing  to  obtain  reliable  data  on  the  material's 
properties  of  large  subsamples.  Lastly,  designers  must  carry  out 
analysis  in  as  much  detail  as  possible  so  that  normal  work  can 
be  carried  out  on  structures  even  under  the  worst  conditions. 
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COMPUTATION  OF  STRENGTH  AND  STABILITY  FOR  A  BUILT  STRUCTURE 

by  Cai  Yiping,  Song  Zenghao,  Wang  Zhihua,  Zhang  Xianwu 
(Computation  Technique  Research  Institute  of  the  China  Aviation 
Research  Institute) 

Zhang  Zugian 

(Design  Institute  of  the  Hongan  Aircraft  Company) 

Abstract 

This  paper  is  devoted  to  solving  strength  and  stability 
problems  of  a  built  structure  by  using  the  finite  element 
method  and  the  increment  method  for  the  geometrical  nonlinear 
structural  problem.  Structural  elements  include  bar,  plane, 
beam,  space  beam,  triangular  and  rectangular  plates  with  mem¬ 
brane  stresses,  rectangular  plate  element  with  bending  deforma¬ 
tions.  Harmonization  between  platebar  and  beam  is  carried  out 
by  using  infinitely  stiff  component.  When  the  generalized  eigen¬ 
value  problem  is  solved  with  the  aid  of  a  simultaneous  itera¬ 
tion  algorithm,  the  initial  test  vector  may  be  automatically 
generated  by  means  of  random  numbers.  In  the  paper  some  tech¬ 
niques  in  computational  procedure  are  also  described,  and  a  few 
numerical  examples  are  given. 

We  have  written  a  general  purpose  program  in  FORTRAN  IV 
language.  The  program  possesses  the  advantages  of  short  comput 
ing  time,  good  generality  and  convenience  for  use. 

I.  Preface 

Among  structural  problems,  when  deflection  is  large  enough  to 
cause  noticeable  changes  in  the  structure's  geometrical 
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configuration  to  the  point  where  an  equilibrium  equation  must 
be  established  based  on  the  configuration  after  deformation, 
geometrically  nonlinear  phenomena  occur.  An  important  condition 
of  geometrically  nonlinear  problems  is  large  displacement  and 
small  strain.  The  properties  after  elastic  buckling  of  the 
structure  is  one  example  of  this  subcategory. 

The  computation  of  the  stability  of  complex  structures  is  an 
important  problem  for  structural  analysis.  Because  of  the  exist¬ 
ence  of  geometric  nonlinearity  the  stress  conditions  of  the  beams 
and  plates  are  generally  Relatively  complex  and  the  boundary  con¬ 
ditions  are  generally  very  difficult  to  determine.  Thus,  to 
attempt  to  find  a  solution  beginning  from  differential  equations 
is  quite  difficult. 

Since  the  1960 's,  the  development  of  the  finite  element  method 
and  numerical  value  analysis  method  have  been  brought  forth  as 
powerful  tools  for  solving  this  problem.  Used  for  linear  struct¬ 
ures,  it  combines  the  developed  matrix  analysis  method  with  the 
gradually  linearized  increment  method,  being  a  means  with  which 
to  solve  nonlinear  structure  analysis  problems.  At  present,  there 
are  many  methods  used  for  solving  the  large  deflection  problems 
of  geometrical  nonlinear  analysis  and  generally  speaking  they  can 
be  divided  into  three  categories:  the  increment  method, itera¬ 
tion  method  and  combination  method.  The  good  and  bad  points  of 
the  methods  have  been  discussed  by  many  writers  f  1]  and  is  not 
within  the  scope  of  this  paper.  Therefore,  we  will  use  the  incre¬ 
ment  method  because  that  method  is  relatively  intuitive,  is  con¬ 
venient  to  use  and  is  able  to  entirely  describe  the  load- 
deflection  characteristics. 

Mathematically,  computation  of  stability  requires  solving  a 
generalized  eigenvalue  problem 

KEq+  A  KQq=0 
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In  the  formula,  K„  is  the  total  elastic  stiffness  matrix  and  K_ 

Cj  G 

is  the  total  geometric  stiffness  matrix.  Based  on  the  model's 
smallest  eigenvalue  we  can  compute  the  structure's  destabilized 
critical  load.  Further,  the  corresponding  characteristic  vector 
is  a  destabilized  waveform. 

We  used  FORTRAN  IV  language  to  write  a  general  purpose  program 

for  computing  the  strength  and  stability  of  the  plate,  bar  and 

beam  built  structure  as  well  as  to  compute  some  examples  and 

problems.  The  results  were  relatively  satisfactory. 

/ 

II.  Computation  Methods 

1.  When  an  elastic  continuous  medium  is  in  a  small  deflection, 
the  strain-displacement  equation  can  be  viewed  as  linear  and 
when  existing  within  external  force  p  and  displacement  U,  there 
is  the  following  linear  relation: 

KU-p  (1) 

Yet,  when  the  structure  has  large  deflection,  we  must  consider 
the  geometric  configuration  changes  caused  by  displacement.  At 
this  time,  the  high  order  item  (nonlinear  item)  in  the  strain- 
displacement  equation  cannot  be  overlooked.  These  nonlinear 
items  cause  the  element's  stiffness  matrix  k  to  be  changed  into 

*-*«  +  *■  (2) 

K_  is  the  ordinary  elastic  stiffness  matrix  and  k_  is  the  geo- 

£j  G 

metric  stiffness  matrix.  They  are  not  only  related  to  the  geo¬ 
metric  measurements  but  also  connected  to  the  internal  force 
existing  when  the  computation  steps  begin.  After  the  elastic 
and  geometric  stiffness  matrices  are  ascertained,  we  can  form 
the  total  stiffness  matrix 

K-Kt+K9  (3) 
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In  this  way,  we  can  use  a  conventional  method  to  solve  the  dis¬ 
placement  and  stress. 


2.  To  consider  the  geometric  configuration  changes  when  the 
external  load  gradually  increases,  we  used  a  series  of  linear 
steps  to  make  approximations  for  the  nonlinear  problem.  Each 
step  represented  one  load  increase.  This  method  is  called  the 
increment  method  £1,2]. 


While  using  the  increment  method  procedure,  each  step  should 
compute  the  model’s  nodal  point  coordinates  based  on  the  dis¬ 
placement  increased  changed  structure  and  then  proceed  to  find 
the  new  stiffness  matrix  and  geometric  matrix. 

3.  When  a  certain  metric  standard  p*  is  selected  for  the  ex¬ 
ternal  load,  the  external  load  is 

★ 

P=  A.  p 


Because  the  geometric  stiffness  matrix  forms  a  ratio  with 
the  internal  force  when  the  steps  begin,  thus 


Therefore  we  have 


U  x/^sr'ip* 


We  can  see  from  this  that  when 


l/C.+XKSI-O 


(4) 


displacement  U  tends  to  be  infinitely  great.  We  found  the  modular 
minimum  value  of  A  from  equation  (4)  and  represented  it  as 


A 


critical ’ 


Thus  we  can  obtain  the  structure's  destabilized 
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critical  load 


P%»  “!«»■/**  (5) 

Key:  1.  Critical?  2.  Critical 
It  can  be  seen  from  equation  (4)  that  the  structure's 
destabilization  equation  is 

(K,+\Ki)  g  =.  o  (6) 


This  is  a  generalized  eigenvalue  problem.  The  minimum  eigen¬ 
value  found  is 

vector  is  a  destabilized  waveform. 


^“critical  an<*  ^ts  corresPon<^in9  characteristics 


4.  The  structural  elements  include  the  bar  triangular  and 
rectangular  plates  rectangular  plate  element  with  bending  deform¬ 
ations,  beam  and  space  beam.  Their  elastic  stiffness  matrix  and 
geometric  stiffness  matrix  can  be  seen  in  references  [3,4,5].  To 
economize  on  the  contents  and  extend  the  scale  of  computations, 
all  of  the  elastic  stiffness  matrices  and  geometric  stiffness 
matrices  after  total  loading  used  a  one-dimensional  change  band 
width  store  pattern. 

5.  The  harmonization  between  the  beam,  plate  and  bar  used  the 
infinite  stiffness  element  process  method.  The  stiffness  compon¬ 
ent  used  in  the  joining  area  of  the  plate,  bar  and  beam  can 
guarantee  that  the  deformations  in  each  element  are  harmonized. 
The  use  of  the  rigid  body  displacement  law  can  determine  the 
displacement  relation  of  each  nodal  point  on  the  infinite  stiff¬ 
ness  component  and  the  use  of  the  principle  of  virtual  work  can 
find  the  relation  of  each  nodal  point  on  the  rigid  body.  Here, 

we  will  only  give  as  an  example  the  joining  of  bar  jk  and  beam 
as  showing  in  fig.  1.  Other  cases  can  be  attained  analogously 
[5]  . 
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j£#ff  (2) 


XRfMttScO) 


(4) 

(5) 


Fig.  1  A  Case  of  Joining  Bar  to  Beam 

Key:  1.  Basic  standard  coordinate  system 

2.  Element  bar 

3.  Infinite  stiffness  element 

4.  Beam  axis  line 

5.  Local  coordinate  system 

Here,  component  ji  is  an  infinite  stiffness  component, 
point  i  is  the  fixed  joining  point,  point  j  is  the  hinge  point, 
and  the  displacement  of  joint  j  is  determined  by  the  dis¬ 
placement  and  corner  of  point  i,  that  is 


In  the  formula 


r  o 


—  C03  Y 


*■>  cop  3 


COSY  -cos  3  ' 
0  cos  a 
—coso  o 
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is  the  third  order  unit  matrix,  is  the  length  of  stiff¬ 

ness  component  ji,  and  cos  a  ,  cos  0  and  cos  Y  are  the  direction 
cosines  of  ij  . 


If  the  jk  bar's  original  stiffness  matrix  is  after 

using  point  i  displacement  to  replace  point  j  displacement,  the 
stiffness  matrix  of  bar  jk  is 


Tst 


T  .  [Tit  ' 
I  Cku)$M§l 


(8) 


The  equilibrium  equation  is 


f0" 

f  u,' 

Qw 

Vi 

0./ 

Wl 

M.I 

9., 

M « 

1  “ 

0„ 

M.i 

9., 

Q* 

Q,» 

■  Qmt 

This  type  of  point  j  does  not  appear  in  the  total  stiffness 
matrix  and  in  the  same  way  the  external  load  on  point  j  must 
also  have  an  equivalent  transformation  up  to  point  i 

<P,>  -CTuVip,}  (lo) 

6.  To  obtain  the  characteristic  value  and  characteristic 
vector  as  well  as  the  partial  characteristic  value  and  charact¬ 
eristic  vector  of  equation  (6) ,  we  used  a  simultaneous  iteration 
algorithm  [7],  At  this  time,  equation  (6)  can  be  written  in  the 
following  forms: 
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or 


Bx-HAx  (11) 

★ 

In  the  formula,  A(e.a.  K_)  and  B(e.g.-K_)  are  both  real  symmetrical 

matrices  and  A  is  fixed.  Using  Cholesky  decomposition,  A=LLT, 
equation  (11)  then  changes  into  the  ordinary  characteristic 
value  problem  of  the  real  symmetrical  matrix 


.  LrlBL~ry-*y 

(12) 

T 

In  the  equation  y=L  x 

(13) 

Several  characteristic  values  of  are  the  largest  derived 

by  a  simultaneous  iteration  algorithm  but  we  required  the  small¬ 
est  1'.  Thus,  when  using  this  method  we  exchanged  the  K_  and 
*  k 

Kq  positions. 

It  should  be  pointed  out  that  before  using  a  simultaneous 

iteration  algorithm,  we  must  first  carry  out  matrix  compression 
* 

of  K_  and  K_,  and  delineate  the  lines  and  rows  corresponding 
E  G 

to  the  known  displacements.  This  is  more  complex  than  using  the 

★ 

one-dimensional  change  band  width  stored  K„  and  K„.  If  the 

E  G 

method  is  not  suitable  and  the  order  of  the  matrix  is  high,  ex¬ 
penditures  will  be  startling.  Our  designed  "revised  fixed 
position  method"  computed  and  determined  the  positions  of  each 
element  after  delineating  the  lines  and  rows  from  back  to  front 
so  that  it  was  only  necessary  for  each  element  to  shift  once. 

This  greatly  economized  on  time. 
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It  is  worth  pointing  out  that  if  a  person  hopes  to  find  X 

characteristic  values  then  it  is  much  better  to  use  m  (  >1  > 

vectors  for  iteration.  This  is  because  the  converging  speed  of 

k  characteristic  vectors  depend  on  iu,  ,  and  the  smaller 

m  k 

the  specific  value  the  faster  the  convergence.  Taking  m  ^  X  > 
the  extra  amount  of  work  expended  for  each  time  of  iteration  is 
compensated  for  because  of  the  large  raising  of  convergence 
speed.  Experience  has  proven  that  taking  m  as  about  1.5  times  of 
X  is  relatively  suitable. 

Finally,  we  think  that  a  simultaneous  iteration  algorithm 
not  only  improves  convergence  but  also  does  not  necessitate  any¬ 
thing  like  the  power  method  wherein  after  each  time  of  finding  a 
characteristic  value  it  is  necessary  to  carry  out  contraction  of 
the  matrix.  When  an  initial  approximation  is  provided  for  the 
characteristic  vector,  the  computer  program  still  possesses  the 
function  of  using  random  numbers  to  produce  initial  test  vectors. 

III.  Computation  Examples 

Example  1.  On  p.  268  of  reference  [2]  there  is  an  analysis 
of  the  stability  of  a  simple  truss.  Our  computation  results  are 
the  same  as  those  of  reference  [2] . 

Example  2.  We  also  carried  out  stability  analysis  of  a  column 
v.ith  one  end  solidly  joined  and  the  other  end  hinged  which  was 
found  on  p.  271  of  reference  [2] .  Computation  results  were 
identical . 

Example  3.  We  analyzed  the  stability  of  the  plate  and  bar 
built  structure  shown  in  fig.  2.  The  structure  is  composed  of  20 
bars,  6  triangular  plates  and  8  rectangular  plates.  There  are  a 
total  of  12  nodal  points:  the  4  nodal  points  on  the  lower  part 
are  fixed  and  the  4  on  the  top  part  separately  sustain  external 
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loads  of  1,000  kilograms.  The  bar  element's  sectional  area  is 
2 

A=1  cm  and  the  plate  element's  thickness  is  t=0.1  cm.  Elastic 

6  2 

modulus  E=0. 72x10  lg/cm  ,  Poisson’s  ratio  >U»=0. 3,  and  the  load 

increments  are  0.1*  0-3*,  0.3*,  0.3*,  respectively. 

P  P  p  P 


Fig.  2  A  Structure  Built  of  Plates  and  Bars 


After  computations  we  obtained 

A-  . _  •  .=115.3 

''critical 

Thus,  on  the  top  four  points  there  is 


pcritical=115,3xl°3  kil°grams 


Computation  results  prove  that  the  distribution  of  displace¬ 


ment  and  stress  conforms  to  the  pattern,  and  the  A. 
value  and  related  materials  match. 


critical 


IV.  Conclusion 


1.  In  essence,  the  stability  problem  of  a  built  structure 


should  belong  in  the  nonlinear  theoretical  range  of  elastic 
mechanics.  In  computations,  we  used  the  microdisturbance  criteria 
of  the  nonlinear  large  deflection  theory  and  elastic  stability 
theory.  A  built  structure  is  not  a  perfect  system  and  a  large 
quantity  of  computations  have  shown  that  the  critical  loads  of 
built  structures  are  extreme,  when  a  built  structure  reaches 
critical  load,  the  level  of  the  structure's  main  support  force 
component  stress  is  relatively  low  but  the  levels  of  stress  of 
certain  secondary  components  is  relatively  high.  Thus,  proper 
adjustment  of  the  geometric  measurements  of  components  with 
higher  stress  levels  is  advantageous  to  raising  the  structure’s 
critical  load. 

2.  When  using  the  increment  method  for  solution,  the  external 
loads  make  up  gradual  step  loading.  During  each  step  nonlinear 
influence  is  small  and  it  is  only  necessary  that  the  load  in¬ 
crement  be  sufficient.  The  solution  of  each  step's  increment  can 
be  obtained  accurately  and  the  precision  of  the  displacement's 
total  value  and  internal  force  can  be  raised.  In  reality,  even 

if  we  use  the  relatively  large  linear  computation  steps  of  these 
problems,  we  will  still  obtain  results  very  close  to  the 
structure's  real  nonlinear  properties. 

Because  the  first  step  of  the  increment  method  is  an  elastic 
solution,  the  added  load  increment  must  be  as  small  as  possible 
so  as  to  decrease  the  computational  errors  of  the  displacement. 

3.  After  testing  several  examples,  computations  were  carried 
out  for  a  certain  fuselage  composite  beam  computation  model  and 
the  results  were  satisfactory. 
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STUDY  ON  TIME  BETWEEN  OVERHAULS  (TBO)  OF  AIRCRAFT  TURBOJET 
AND  ITS  RELIABILITY 

by  Zhang  Yimin 

(Shanghai  Aviation  Industry  Corporation) 

Abstract 

In  this  paper  statistical  methods  (frequency  histogram, 

2 

probability  paper,  linear  regression,  chi-square  ~K  tests) 
are  applied  to  analyzing  TBO  of  aircraft  engines.  The  calcu¬ 
lated  results  of  three  types  of  aeroengines  have  shown  that  the 
normal  distribution  is  an  appropriate  model  for  TBO  of  aircraft 
engines.  Therefore,  along  with  the  given  value  of  TBO,  it  is 
reasonable  to  indicate  the  probability  of  the  actual  operating 
time  of  aircraft  engines  exceeding  the  given  TBO.  It  is  sug¬ 
gested  to  express  the  TBO  by  HfR^*).  Thus  many  unnecessary 
overhauls  can  be  avoided  which  is  of  significant  economical  im¬ 
portance.  The  mean  TBO  f  normal  distribution  if=H(50%)  is 
recommended  as  the  parameter  of  aircraft  engines  for  overhaul. 
While  accepting  the  mean  TBO  concept,  on-condition  maintenance 
and  condition  monitoring  technology  must  be  adopted  in  engine 
operation  in  order  to  obtain  sufficient  flight  safety.  If  the 
over  50%  actual  operating  time  of  aircraft  engines  is  greater  than 
the  mean  TBO  in  airline  operations,  the  TBO  can  be  prolonged. 

The  method  for  estimating  the  prolonged  life  is  also  given. 

Reliability  is  obviously  very  important  in  the  operating  per¬ 
formance  of  aircraft  engines.  This  is  because  the  reliability  of 
an  aircraft  engine  is  directly  related  to  the  flight  safety.  At 
the  beginning  of  the  1950' s  when  turbojets  were  first  put  into 
operation,  very  stringent  TBO  control  methods  were  adopted  to 
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guarantee  the  reliability  of  engine  operation.  Yet  operation 
experiences  showed  that  TBO  could  not  be  an  index  for  the 
reliability  of  flight  operation  [1]  and  that  after  a  large  number 
of  engines  operated  in  excess  of  the  TBO  they  could  still  operate 
normally.  However,  many  engines  had  breakdowns  within  the  stipul¬ 
ated  TEO  period  and  moreover  the  reliability  of  the  engines  did 
not  increase  with  operation  time  to  a  certain  value  but  notice¬ 
ably  worsened.  Therefore,  we  indicated  the  reliability  with  prob¬ 
ability,  within  a  fixed  amount  of  time  and  under  stipulated 
operation  conditions,  there  was  no  hindrance  to  giving  play  to 
the  probability  ofp  stipulated  performance.  Fig.  1  gives  the  fre¬ 
quency  histogram  of  the  TBO  distribution  for  284  CF6  engines 
within  917,500  flight  hours.  The  longest  reached  to  over  6,000 
hours  and  the  shortest  was  only  100  odd  hours.  The  TBO  mean  value 
H  was  3,231  hours  and  variance  s  was  1,358  hours.  Fig.  1  also 
qives  the  curve  calculated  according  to  the  normal  distribution, 
Weibull  distribution  and  exponential  distribution  which  shows 
that  the  CF6 ' s  TBO  distribution  is  close  to  the  normal  or  Weibull 
distributions . 


Fig.  1 


91 


F/G  12/1 


AD-A118  973  FOREIGN  TECHNOLOGY  OXV  WRIGHT-PATTERSON  AFB  OH 
ACTA  AERONAUTICA  ET  ASTRONAUTICA  SlNICA.tU) 

JUL  82 

UNCLASSIFIED  FTb-ID<RS> T-0518-82  NL 


Fig.  1  Frequency  Histogram  For  TBO  Data  of  CF6 

Key:  1.  Line  1  -  normal  distribution 

2.  Line  2  -  Weibull  distribution 

3.  Line  3  -  exponential  distribution 

4.  Number  of  overhauled  engines  M. 

5.  TBO  H  (hours)  K 


Further,  the  operation  reliability  and  TBO  relation  of  the 
CF6  is  separately  drawn  on  normal  distribution  probability  paper 
logarithmic  coordinate  paper  and  Weibull  distribution  coordinate 
paper.  Operation  reliability  is  the  percentage  occupied  in 
a  normally  operating  engine  within  stipulated  time  Hv.  It  is 
also  the  engine  percentage  when  the  TBO  exceeds  hours.  There¬ 
fore  : 


/?*=*!  — 


■ » i 


N 


(1) 


In  the  formula,  M  is  the  TBO's  engine  number  in  the  i  interval 
and  N  is  the  total  number  of  operating  engines  within  this  batch 

We  used  linear  regression  analysis  for  the  24  groups  of 
statistical  data  for  the  CF6  in  fig.  1  to  find  its  regression 
function  formula.  In  this  way,  there  were  many  errors  yet  calcul 
at ions  were  convenient,  could  be  carried  out  on  a  programmed 
calculator  (such  as  TI-59)  and  need  not  use  a  computer.  For 
standard  normal  distribution  N  (0,1),  standard  deviation  Z  was 
the  variable  [2] ,  and  the  definition  of  Z  is: 


In  the  formula,  &  is  the  generating  variance.  Yet,  when  the  TBO 
is  based  on  the  normal  distribution,  the  relation  between  relia¬ 
bility  and  Z  is 
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*«-  1  -4>(Z,)  = 


_ 1_ 

✓  2  * 


Z*_ 

<(Z 


(3) 


In  the  formula,  ^  (Z^)  is  called  the  normal  probability  integral 
and  we  can  table  look-up  [2]  or  use  a  program  stored  in  the 
TI-59  calculator  for  calculations. 

The  23  groups  of  data  in  CF6  (see  fig.  1,  the  one  group  of 
R^=0  is  deleted)  finds  R  based  on  equation  (1)  and  then  from  the 
derived  we  find  by  using  eauation  (3) .  We  then  derive  23 
groups  of  data  for  .fT  and  Z  and  by  regression  analysis  obtain: 

Z --2. 43175  +  0. 00077158//  (  4  ) 

r  -  +  0.9876 

The  correlation  coefficient  is  r=+0.9876.  Equation  (4)  and  the 
23  groups  of  data  (H^ , Z^)  are  drawn  on  normal  probability  coordin¬ 
ate  paper.  We  know  from  fig.  2  that  the  23  groups  of  data 
(H^ , Z^)  are  very  close  to  regression  line  (4)  and  that  when  the 
correlation  coefficient  is  close  to  1  the  TBO  of  CF6  coincides 
with  the  pattern  of  normal  distribution. 
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Fig.  2  Estimated  Regression  Line  According  to  Normal  Distribution 
and  Reliability  of  CF6 

Key:  1.  Reliability  R^  (%) 

2.  TBO  H  (hours) 

3.  Equation  (2) 

When  we  used  the  same  regression  analysis  method  (23  points) 
for  the  exponential  distribution,  we  obtained: 

!g/?  -2. 44308-0. 00030828H 

When  correlation  coefficient  r=  -0,83821,  this  shows  that  the 
linear  relation  between  lgR  and  H  is  not  close  and  from  fig.  3 
we  know  that  the  actual  measuring  points  are  far  apart  from  the 
regression  line. 
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Fig.  3  Estimated  Regression  Line  According  to  Exponential 
Distribution  and  Reliability  of  CF6 

Key:  1.  Reliability  R4  (5) 

2.  TBO  H  (hours) 

3.  Equation  (3) 

For  the  Weibull  distribution,  we  used  regression  analysis  to 
determine  configuration  parameter  m  and  measurement  parameter  Sq 
we  determined  position  parameter  aQ  based  on  the  in  the 


A  test  which  is  the  minimum  value  condition, 
mined  by  the  following  equation: 


is  deter - 


(NPt-M>y 

NP> 


In  the  equation,  is  the  probability  of  the  operating  TBO  in 
the  k  interval.  See  table  1  for  the  calculated  results  of  the 
A  tests.  It  was  more  suitable  for  to  be  taken  as  -1200. 


| 

I 

(4  J  *  *  * 

»  * 

n*  ft 

i 

1  <H  '  0 

o»  *  - 1100 

o»  «  -  1200 

j  at  •  -  130* 

<5)»**Scr 

24 

0.99537 

0.99531 

0.9*534 

x*-  2  x* 

k  -  1 

32.247  | 

32.24* 

32.  It* 

Table  1  The  Calculated  Results  of  /2  Test 

Key:  1.  Type  of  distribution 

2.  Normal  distribution 

3.  Exponential  distribution 

4.  Weibull  distribution 

5.  Correlation  coefficient  r 


In  this  way,  we  find  that  the  Weibull  distribution  regression 
formula  is: 


In  ln-^' -  -  30 . 9575  +  3 . 647861n(/f  4- 1200)  (  7  ) 

When  correlation  coefficient  r=0. 99538,  this  shows  that  its 
linear  degree  of  correlation  is  relatively  good.  We  know  from 

fig.  4  that  each  actual  measuring  point  is  relatively  close  to 
the  line  (equation  7) . 
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Fig.  4  Estimated  Regression  Line  According  to  weibull 
Distribution  and  Reliability  of  CF6 

Key:  1.  (Hours) 

2.  Equation  (7) 

Therefore,  we  know  from  regression  analysis  that  the  real  dis¬ 
tribution  and  normal  distribution  or  Weibull  distribution  of  CF6 
are  all  quite  close. 

To  carry  out  the  X  test  for  the  normal  distribution,  we 
took  samples  of  H  and  s  as  the  and  Cf  of  the  generating  pop¬ 
ulation,  took  significance  a  as  0.05  (e.g.  believability  was 
0.95)  and  degree  of  freedom  v=n-3=21  so  *".  to  satisfy  the  follow¬ 
ing  conditions  [2] : 


Xj._«_<X*<Xi,  1  (8) 
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Then,  the  original  supposition  could  be  accepted  and 

y  2  _in  to,  y  2  =35.479.  We  know  from  table  1 

/'  2 1.0.025“1U*  ''21 ,0.975 

that  both  the  normal  distribution  and  Weibull  distribution 
satisfy  the  conditons  of  (8)  and  this  can  accept  the  supposition 
of  the  normal  distribution  and  Weibull  distribution.  When 
X  e=262.81  during  exponential  distribution,  this  does  not  sat¬ 
isfy  the  conditions  of  (8)  and  thus  eliminates  the  supposition 
of  the  exponential  distribution. 

Further,  taking  the  JT3D-7  engine  as  an  example,  there  was  a 
total  of  224,798  flight  hours  from  July,  1974  to  August,  1980. 

The  greatest  number  of  flight  hours  for  one  aircraft  was  8,039 
hours  while  the  manufacturing  plant's  recommended  TBO  was  9,000 
hours.  According  to  the  old  TBO  concept,  it  is  not  necessary  to 
dismantle  an  engine  for  overhaul.  Yet  during  this  period  of 
actual  operation  4  engines  stopped  in  midair  and  34  were  changed 
early.  The  midair  stoppage  rate  was  0.0178  and  the  early  change 
rate  was  0.165.  The  total  flight  time  of  these  38  engines 
between  overhauls  was  127,221  hours;  there  were  4  which  did  not 
operate  up  to  1,000  hours  before  having  to  be  dismantled  for 
overhaul  and  the  longest  time  was  7,296  hours.  The  38  engines 
were  divided  into  7  intervals,  H  was  3,  316  hours  and  variance 
was  1,798  hours.  According  to  the 

X  2  test  carried  out  for  normal  distribution,  X  2=  A 

.2  2  *«i 

Xy.=5.771.  We  know  from  reference  [21  that  X4  g  q25=^'484  an(* 

2  2 

X  4  Q  025=0.484  and  ^4  0975=11,143  which  satisfy  the  conditions 

of  (8)  and  therefore  can  accept  the  supposition  of  normal  distri- 

2 

but ion.  Based  on  X  as  the  minimum  value  condition,  we  found 
the  Weibull  distribution  regression  formula; 
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Inin  ^  =-22. 9449  +  2. 68287la(//*+ 1300)  (  9  ) 

2  <2  2 

Correlation  coefficient  r=0. 99345  and  pf  =i<  satisfy 

K*  I  K 

the  conditions  of  (8)  and  can  therefore  also  accept  the  supposi¬ 
tion  of  the  Weibull  distribution. 

Statistics  on  55  TBO  for  AH-20M  have  shown  that  because  of 
the  influence  of  the  old  TBO  concepts,  after  the  operation  of 
19  of these  engines  reached  the  stipulated  TBO  of  3,000  hours 
they  were  dismantled  and  sent  for  overhaul  without  having  broken 
down.  Thus,  testing  should  be  carried  out  using  one-sided 
truncated  normal  distribution  and  linear  regression  analysis 
should  be  used  to  determine  the  mean  TBO  and  variance: 

Z --1. 84160+0. 000768046/f.  (10) 

The  correlation  coefficient  is  r=0. 99581.  By  comparing  formulas 

(2)  and  (10),  we  can  know  that  i  =1/0.000768046  and  H  =1.84169. 

°  a 

Therefore,  when  the  variance  is  1,302  hours,  the  mean  TBO  is 

7 

2,398  hours  and  we  calculate  X^=  2 

k 

which  satisfies  the  conditions  of  (8) .  Thus  it  can  accept  the  one¬ 
sided  truncated  normal  distribution  supposition.  If  the  truncated 
supposition  is  not  used,  even  though  we  consider  the  TBO  of  the 
19  dismantled  3,000  hours,  H  os  2.136  hours,  s  is  890  hours,  and 
2  2 

c  z  Pf  ^=19.208  which  do  not  satisfy  the  conditions  of  (8). 

Thus,  it  is  necessary  to  use  the  one-sided  truncated  normal  dis¬ 
tribution.  In  the  same  way,  we  also  use  one-sided  truncated 
Weibull  distribution  for  testing.  Using  regression  analysis  we 
can  obtain: 
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In  In— 5—  -  -  23 . 0046  +  2 . 783831n  (N,+ 1 000 )  (11) 

2  7  2 

Correlation  coefficient  r-0. 99538  and  ft  =  Xv=3.617. 

k=l  K 

Thus,  it  can  accept  the  one-sided  truncated  Weibull  distribution. 
The  above  calculations  can  be  completed  on  TI-59  programmed 
calculator  programming  and  a  stored  calling  program. 

Based  on  the  TBO  statistical  tests  of  the  above  three  types  of 
engines,  the  normal  distribution  and  Weibull  distribution  are 
both  acceptable.  Yet  m  is  about  3  in  both  the  normal  and  Weibull 
distribution  so  that  the  use  of  normal  distribution  was  conven¬ 
ient  for  analysis  and  calculations  and  therefore  normal  distribu¬ 
tion  was  used  in  later  analysis.  It  can  be  seen  from  the  above 
examples  that  changes  in  TBO  concepts  have  caused  great  changes 
in  engine  operation  and  maintenance. 

1.  According  to  the  old  TBO  concept  of  the  1950's,  even  if 
the  stipulated  TBO  for  the  JT3D-7  was  1,000  hours,  10.5%  of  the 
engines  would  not  reach  the  stipulated  TBO  but  would  have  to  be 
dismantled  and  overhauled  by  the  manufacturing  plant.  Because 
the  other  34  engines  would  have  already  reached  TBO  by  1,000 
hours  of  operation,  they  also  had  to  be  dismantled  and  overhauled. 
This  caused  the  34  engines  to  have  92,388  fewer  flight  hours 
occupying  73.1%  of  the  total  flight  time  of  126,388  hours  for  the 
34  engines.  This  is  unnecessary  waste  and  shows  that  the  use  of 
the  new  concepts  possesses  significant  economic  worth. 

2.  Along  with  providing  the  TBO  we  should  also  give  the 
probability  of  the  operating  life  reaching  the  indicated  TBO. 

The  higher  the  requirements  for  reliability  the  shorter  the  TBO 
(see  table  2) . 
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Ci)  - - 

JT3D-7 

CFO 

AH-zoM 

(3)  F.-sXttna***.  Hi 5> 

6300 

S2S0 

>3000 

(4)fl*-50 **«■«**.  Hi so> 

SOM  ; 

SIM  ' 

2300 

C5)  Hi 90) 

•  1000  1 

1150 

625 

Table  2  The  Relationship  Between  TBO  and  Reliability  For 
Several  Aeroengines 

Key:  1.  TBO 

2 .  Engine  Model 

3.  Rk=5%  of  the  time  TBO,  H(5) 

4.  Rk=50%  of  the  time  TBO,  H(50) 

5.  Rk=90X  of  the  time  TBO,  H(90) 

Therefore,  the  use  of  H(Rk%)  is  proposed  for  indicating  TBO  and 
we  can  use  the  mean  TBO  to  replace  the  old  TBO.  However,  if  the 
relationship  between  the  probability  of  the  actual  operating  life 
reaching  the  mean  TBO  and  its  distribution  is  subordinated  to 
normal  distribution,  its  probability  is  50%  and  thus  one  half 
of  the  engines  can  reach  the  mean  TBO.  Table  3  gives  the  per¬ 
centages  of  the  operating  lives  of  several  engines  which  exceed 
the  mean  TBO.  All  approach  50%.  Table  3  also  gives  TBO  statistics 
for  several  engines  used  by  the  United  States  military  [31  which 

shows  that  the  above  new  concept  for  TBO  is  also  applicable  in 
engines  for  military  use. 
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(1)  (**-«*) 

CFe 

(DC-10) 

JT3D-7 
(Ktf  707) 

AH-24 

AH-20M 

J57-P19 

(B-52D) 

J57-P21 

(F-100) 

TFS0-P-3 

(F-lll) 

TF41-A1 

CA-7D) 

284 

102 

52 

113 

187 

137 

115 

¥***.  (4) 

1231 

3348 

3099 

2396 

3666 

721 

526 

333 

(5 

)  41.62 

42.11 

51.96 

47.27 

)  6122 

7286 

>4000 

>3000 

4000 

1000 

1000 

750 

Table  3  The  Mean  TBO  of  Several  Aeroengines 

Key:  1.  Engine  model 

2.  (Boeing  707) 

3.  Statistical  numbers 

4.  Mean  life,  hours 

5.  Exceeding  mean  life, 

6.  Highest  TBO,  hours 

3 .  The  use  of  the  midair  stoppage  rate  as  an  aircraft  re¬ 
liability  index.  After  using  the  mean  TBO  concept,  about  one 
half  of  the  aircraft  still  did  not  reach  the  mean  TBO  but  had 
breakdowns.  Therefore,  within  each  1,000  engine  flight  hours 
the  number  of  times  the  engine  itself  had  midair  stoppage  and 
the  number  of  times  it  had  to  be  changed  early  are  called  the 
midair  stoppage  rate  and  early  change  rate.  We  used  the  midair 
stoppage  rate  as  the  main  engine  reliability  index.  Engine 
manufacturing  plants  carry  out  overhauling  of  engines  for  early 
changes  but  are  not  responsible  for  inclusive  repairs.  They  are 
only  responsible  for  the  proportional  compensation  of  important 
parts  damage  owing  to  design  and  manufacturing  during  the 
guarantee  period  but  the  parts  guarantee  period  is  obviously 
less  than  their  life.  For  example,  the  operation  limit  of  a 
JT3D  compressor  disc  can  reach  20,000  hours  yet  the  guarantee 
period  is  only  4,000  hours. 
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4 .  On-condition  maintenance  and  condition  monitoring  must 
be  used.  The  use  of  the  mean  TBO  concept  signifies  that  one 
half  of  the  engines  will  reach  the  aforementioned  requirement 
changes  of  the  mean  TBO.  Thus,  if  there  is  a  breakdown  and  the 
engine  does  not  reach  the  designated  time  for  changes,  there 
will  be  midair  stoppage  or  delayed  error  flight.  The  approx¬ 
imate  relation  can  be  indicated  by  the  following  formula: 

IFSD~5™-PEX  (12) 

H 

In  the  formula,  IFSD  and  PER  separately  represent  the  midair  , 
stoppage  rate  and  early  change  rate.  Thus  it  is  necessary  to  use 
cn-condition  maintenance  and  the  condition  monitoring  technique 
which  is  not  dismantling  and  overhauling  engines  according  to 
a  stipulated  time  limit  but  monitoring  the  engine's  important 
parameters  as  well  as  adopting  appropriate  measures  to  monitor 
the  engine's  key  parts  under  undecomposed  engine  conditions  and 
make  on-condition  changes  of  parts  when  necessary.  As  regards  the 
stipulated  total  life  and  total  rotational  life  of  discs  and 
axes,  these  techniques  can  promptly  discover  possible  breakdowns 
and  utilize  early  change  to  guarantee  flight  safety. 

5.  The  expectation  of  the  mean  TBO.  According  to  che  statis¬ 
tical  materials  on  navigation  operations,  over  50%  of  the 
engines  exceed  the  stipulated  TBO.  The  mean  TBO  can  then  be 
estimated  by  the  following  formulas: 


Z | // .  (  A’. )  —  Z.  / ( , ( /V , ) 

z,-z. 


(13) 


o 


z,-4 


(u) 


In  the  formulas,  Z ^  and  Z 2  are  based  on  the  and  R2  numerical 
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•r~  *  *  nrnr»wi  1  r  r  -r  ■  ■ 


values  determined  from  formula  (3) . 


The  calculation  results  from  formulas  (13)  and  (14)  are 
listed  in  table  4. 


( i )  Hifrtt 

1 

H ( 2 

) 

!  *• 

Zt 

1 

1000 

04. 1 

WEM 

30.85 

-  1.0 

+  0.5 

1007 

667 

1000 

84. 1 

38fl 

St. 2 

-  1.0 

+  0.5 

1700 

700 

1000  1 

07. 7  ! 

MEM 

04.1 

-2.0 

-1.0 

3000 

1000 

Table  4  The  Estimated  Mean  TBO  According  to  Operating 
Statistical  Information  of  Aeroengines 

Key:  1.  hours 

2 .  H2  hours 

3 .  H  hour  s 

4 .  &  hours 

To  decrease  estimation  errors,  the  statistical  engine  numbers 
should  be  larger  than  100  and  generally  should  be  greater 
than  0.5.  At  this  time,  the  error  variance  of  test  statistical 
value  Rj  is  smaller  than  one-tenth  of  Rj .  If  there  is  statistical 
data  for  more  than  two  points,  we  can  then  use  the  linear 
regression  analysis  method  to  estimate  H  and  d . 

Aeroengines  can  also  use  the  prolonged  operating  life  method. 
For  example,  the  United  States  Air  Force  operates  twelve  F100 
engines.  Their  TBO  is  500  hours,  the  continuous  flight  of  two 
F100  engines  exceeded  500  hours  and  TBO  was  extended  to  750  hours 
[4] .  Among  the  operational  N  engines,  the  operating  life  is 
lower  than  the  TBO  of  the  D.  Among  the  spot  checked  n  engines, 
the  probability  of  m  number  of  unqualified  engines  can  be  derived 
by  the  following  formula: 
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p 


(15) 


C(£)C  (?:£)_ 
C(i) 


Table  5  gives  the  probabilities  of  continually  having  two 
engines  exceeding  500  hours  and  based  on  the  maximum  estimated 
D-0~3  they  can  continuously  have  two  engines  exceeding  500 
hours . 


(1)  D 

• 

1 

1  2 

3 

4 

8 

* 

(2)  ******  J?jK 

1M 

91. #7 

13.33 

73 

34.37 

51.33 

59 

(3)  ******  SM*NN«*.  x 

100 

IS. 33 

33.11 

34.3* 

42.42 

31.32 

'22.73 

(4)  -*■*  (M**»  X 

• 

13.17 

39.33 

40.91 

49.43 

53.03 

54.59 

(5)  X 

0 

1.32 

4.96 

141 

15.15 

22.73 

Table  5 


The  Probabilities  of  Two  "Pacer  Century"  Engines' 
TBO  Exceeding  the  Mean  TBO  in  12  Operational 
Aeroengines 


Key:  1.  Life  is  smaller  than  the  number  of 
engines  of  500  hours 


2.  Percentage  with  operating  life  greater 
than  500  hours, 

3.  Probability  of  continuously  having  two 
engines  exceed  500  hours,  % 

4.  Probability  of  one  engine  exceeding 
500  hours  and  one  engine  being  less 
than  500  hours,  % 

5.  The  probability  of  both  engines  being 
less  than  500  hours,  % 


Because  point  estimation  error  is  relatively  great,  to  be  safe 
we  take  D=3  and  if  H (75) =500  hours  from  formula  (2)  we  can 
obtain : 
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H(50) -//(/?♦)-£♦* 


(16) 


In  the  formula,  Z^  is  based  on  which  is  derived  from  formula 

(3)  and  therefore  H (50=H (75) =0 . 67450  so  that  the  engine  can 
prolong  the  life  0.6745  <S  hours.  If  we  take  d  to  be  H/2.33 
which  is  214.6  hours,  we  can  prolong  the  life  145  hours.  There¬ 
fore,  using  probability  analysis  we  can  continuously  operate 
several  engines  in  excess  of  the  mean  TBO  and  this  can  prolong 
the  life.  To  guarantee  flight  safety,  we  should  also  carry  out 
tests  on  key  parts  and  test  and  verify  acceleration  simulated 
test  runs  in  the  manufacturing  plant.  We  should  use  proper 
on-condition  maintenance  and  condition  monitoring  measures  in 
operations  maintenance  to  guarantee  operation  reliability. 
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EXPERIMENTAL  RESEARCH  ON  PERFORATED  ACOUSTIC  LINERS  IN 
TURBOJET  ENGINE  AFTERBURNERS 

by  Ou  Xuebin  and  Ni  Guoxiong 
(Shenyang  Liming  Machinery  Company) 

Abstract 

The  reasonable  choice  of  construction  parameters  for 
perforated  acoustic  liners  in  turbojet  engine  afterburners  was 
investigated  from  the  engineer's  point  of  view.  The  perforated 
acoustic  liner  should  have  as  much  as  possible  high  oscillation 
absorptivity  in  a  quite  wide  frequency  range  by  increasing  the 
volume  of  acoustic  resonant  space  and  the  perforated  area  ratio 
of  the  shield.  The  smooth  geometry  of  the  liner  can  make  the 
secondary  flow  uniform  around  the  outer  chamber  shell  and, 
therefore,  avoid  high  temperature  stripes  proceeding  from  the 
gas  separated-f low  vortex  at  the  wave  valley  on  the  perforated 
acoustic  liner.  Consequently,  the  circular  temperature  difference 
and  the  thermal  stresses  of  the  shell  can  be  decreased,  and  the 
wall  temperature  of  the  afterburner  shell  can  be  kept  down  below 
the  allowable  temperature  of  its  material.  As  a  result,  the 
afterburner  can  operate  reliabily  for  a  long  time. 

Table  of  Major  Symbols 

a  ripple  height  (millimeters) 

b  ripple  span  (millimeters) 

d  diameter  of  small  holes  for  oscillation  absorption 

(millimeters) 

f  intrinsic  frequency  of  perforated  acoustic  liner  (hertz) 

f  gas  pulse  frequency  (hertz) 

G^  ratio  of  secondary  flow  and  total  gas  flow  (percent) 
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h  height  of  wave  valley  (millimeters) 

H  flight  altitude  (kilometers) 

L  length  of  perforated  acoustic  liner  (millimeters) 

1  effective  length  (millimeters) 

M  flight  Mach  number 

n  number  of  small  holes  for  oscillation  absorption  (each) 

A  p  pulse  pressure  (kilograms/centimeter^) 

2 

p  mean  pressure  of  gas  (kilograms/centimeter  ) 

r  radius  of  ripple  curvature  (millimeters) 

t  temperature  (°C) 

time  (seconds 

V  resonant  space  (meters^) 

x  length  of  afterburner  shell  (millimeters) 

a  oscillation  absorptivity 

$  thickness  of  perforated  acoustic  liner  (millimeters) 

0  perforated  area  ratio  (percent) 

&  acoustic  resistance  ratio 

^  acoustic  reactance  ratio 

I.  Preface 

Following  the  raising  of  the  afterburning  ratio  for  turbojet 
aircraft,  the  oil-gas  ratio  in  the  afterburner  continually 
increased  thus  causing  the  afterburner  to  produce  a  considerable 
enlargement  of  the  pulse  amplitude  value.  Moreover,  the 
afterburner's  shell  wall  temperature  became  increasingly  higher 
When  designing  according  to  the  acoustic  resonance  principle, 
when  the  perforated  acoustic  liner  is  placed  in  the  main 
combustion  area  it  can  effectively  inhibit  and  absorb  oscillation 
energy.  It  especially  guarantees  smooth  combustion  in  the  greater 
part  of  the  transverse  high  frequency  range.  The  back  part  which 
does  not  have  a  heat  insulation  shield  and  the  afterburner  shell 
are  directly  exposed  to  the  high  temperature  gas  flow,  and  the 
construction  of  the  perforated  acoustic  liner  has  a  noticeable 
influence  on  the  afterburner's  shell  wall  temperature  [1]. 
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We  installed  various  different  types  of  perforated  acoustic 
liners  in  similar  locations  in  the  same  afterburner.  This 
included  comparative  tests  of  the  perforated  acoustic  liner 
oscillation  damping  temperature  drop  effects  on  a  rippleless 
plate  perforated  acoustic  liner  carried  out  on  ground  rig,  high 
altitude  simulated  test  platforms  and  in  flight.  Fig.  1  is  a 
schematic  of  the  afterburner  construction. 


Fig.  1  Schematic  of  Afterburner  Construction 

Key:  1.  Perforated  acoustic  liner 
2.  Transducer 

II.  Construction  and  Oscillation  Absorptivity  of  Perforated 
Acoustic  Liners 

The  construction  parameters  for  different  types  of  perforated 
acoustic  liners  are  listed  in  table  1.  See.  fig.  2  for  the 
construction  schematic. 
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10 

178  I  17 

318 

5 
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1 
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Table  1 
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Table  I  Construction  Parameters  For  Different  Types  of  Liners 

Key :  1 .  Type 

2.  Dimension 

3 .  Parameter 

4.  Millimeters 

5.  Millimeters 

6.  Millimeters 

7.  Millimeters 

8.  Millimeters 

9 .  Each 

10.  Millimeters 

11.  Meters'* 

12.  Millimeters 


C2) isna*  (3)  >a»«» 


Fig,  2  Schematic  of  Perforated  Acoustic  Liners 

Key.:  1.  Type  III  perforated  acoustic  liner 

2.  Type  II  perforated  acoustic  liner 

3.  Type  I  perforated  acoustic  liner 


For  perforated  acoustic  liners  designed  according  to  acoustic 
resonance  principles,  when  gas  pulse  frequency  f  and  intrinsic 
frequency  f  of  acoustic  resonant  space  in  the  combustion  chamber 
both  produce  resonance,  the  acoustic  reactance  ratio  is  zero. 

At  this  time,  the  perforated  acoustic  liner's  absorption  gas 
oscillation  energy  is  maximum,  that  is,  the  liner's  oscillation 
absorptivity  is  maximum.  When  the  gas  pulse  frequency  deviates 
from  the  resonance  point,  oscillation  absorptivity  gradually 
diminishes.  Fig.  3  gives  the  relational  curve  of  oscillation 
absorptivity  a  and  the  gas  pulse  frequency  of  class  I, II  and  III 
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perforated  acoustic  liners. 


Fig.  3  Relational  Curve  of  Oscillation  Absorptivity  a  and 
Gas  Pulse  Frequency  f 

Key:  1.  fxlO^  (hertz) 


The  curve  is  calculated  by  the  following  formula: 


a 


48 

Te  +  i  )*+x‘ 


(i) 


In  the  formula,  8  is  the  acoustic  resistance  ratio  and  x  is  the 
acoustic  reactance  ratio: 


8 


9.37m 

cJacfrr+T79A?j) 


X 


(2) 

(3) 


"f  is  the  intrinsic  frequency  of  the  resanant  space: 


(4) 
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We  carried  out  tests  on  the  influences  of  the  tested 
perforated  area  ratio  on  oscillation  absorptivity  on  a  ground 
rig  using  an  A  type  engine.  When  the  afterburning  oil -gas  ratio 
was  q^.  =0.057  'v'  0.060  and  when  we  used  type  I  perforated 

acoustic  liners,  the  gas  relative  pulse  pressure  in  the  after¬ 
burner  reached  52.6  percent.  Moreover,  when  there  was  whistling 
and  component  damage,  the  perforated  area  ratio  of  the  perforated 
acoustic  liner  increased  from  0.85  percent  to  1.7  percent,  that 
is  it  doubled.  Working  under  completely  identical  conditions,  the 
relative  pulse  pressure  amplitude  value  in  the  afterburner 
dropped  by  one  half  and  the  whistling  disappeared.  The  measuring 
test  results  match  the  pattern  revealed  in  formula  (1) . 

Generally,  the  oscillation  absorptivity  of  the  turbojet 
engines  is  greater  than  the  0.2  frequency  range  and  this  can 
effectively  inhibit  and  eliminate  the  occurence  of  oscillation 
combustion.  We  call  this  frequency  range  the  effective  oscilla¬ 
tion  absorption  frequency  band.  The  tests  discussed  below 
wherein  there  were  the  antivibration  temperature  drop  effects 
of  the  perforated  acoustic  liners  were  carried  out  on  a  B  type 
engine. 


III.  Test  Results  and  Analysis 

1.  Tests  on  Antivibration  Effects  of  Perforated  Acoustic 
Liners 

We  used  a  pressure  transducer  to  measure  the  gas  pulse 
pressure  on  the  shell  wall  behind  the  perforated  acoustic  liner 
(see  fig.  1)  . 

For  the  pulse  pressure  oscillograms  of  various  types  of 
flight  conditions  with  total  afterburning  see  figures  4-9  and 
for  frequency  spectrum  analysis  see  fig.  10. 
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4  Oscillogram  of  Gas  Pulse  Pressure  in  the  Afterburner 
With  Ty^e  I  Liner 

H=0,  M=0,  p=2.25  kilograms/centimeter2,  =  3.8  percent 

Key:  1.  (Seconds) 


5  Oscillogram  of  Gas  Pulse  Pressure  in  the  Afterburner 
With  Type  II  Liner  * 

H=0,  M=0,  p=2.25  kilograms/cent imeter^ ,  — ^  =  2.5  perent 

P 

Key:  1. (Seconds) 


Oscillogram  of  Gas  Pulse  Pressure  In  the  Afterburner 
with  Type  III  Liner  2  A  D 

=5,  M=1 . 28 ,  p.  2.81  kilograms/centimeter  ,  =  3,7  percent 
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Key:  1.  (Seconds) 
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7  Oscillogram  of  Gas  Pulse  Pressure  In  the  Afterburner 
With  Type  III  Lines 

H=13.5,  M=2.18,  p=1.86  kilograms/cent imeter^ ,  =  4.1 

percent  p 

Key:  1.  (Seconds) 
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8  Oscillogram  of  Gas  Pulse  Pressure  In  the  Afterburner 
With  Type  III  Liner 

H=20,  m=1.53,  p=0.65  kilograms/centimeter^ ,  =6.5  percent 

Key:  1.  (Seconds) 
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9  Oscillogram  of  Gas  Pulse  Pressure  In  the  Afterburner 
With  Type  I  Liner 

H=20,  M=1 . 6 ,  p=0.65  kilograms/centimeter^ ,  =6.5  percent 
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Fig.  10  Frequency  Spectrum  of  Gas  Pulse  Pressure  In  the 
Afterburner  With  Type  III  Liner 

H=20,  M=1 . 6 

Key:  1.  (Decibels) 

2.  (hertz) 

2.  Analysis  of  Results  of  Antivibration  Tests 

Perforated  acoustic  liners  with  designs  based  on  the 
acoustic  resonance  principle  can  inhibit  or  eliminate  pressure 
oscillation  in  the  afterburner,  especially  in  the  high  frequency 
oscillation  range.  Test  results  show  that: 

(1)  Gas  pressure  oscillation  in  the  afterburner  changes  with 
flight  conditions,  and  following  the  gas  mean  pressure  drop  in 
the  afterburner  there  is  a  corresponding  increase  of  the  pulse 
pressure  value. 

(2)  The  increases  of  the  liner’s  acoustic  resonant  space  and 
perforated  area  ratio  raised  the  liner's  absorption,  storage 
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and  dissipation  of  oscillation  energy  and  thus  raised  the 
oscillation  absorption  effects  of  the  perforated  acoustic 
liner. 

(3)  The  perforated  acoustic  linear  used  a  rippled  plate 
construction  which  caused  obstruction  of  pressure  wave 
propagation  and  destroyed  the  resonant  effects  in  the  combustion 
chamber.  Thus,  there  was  an  antivibration  effect.  The  circular 
ripple  construction  helped  to  inhibit  transverse  oscillation 
yet  the  resonant  space  and  perforated  area  ratio  were  close. 

Among  the  perforated  acoustic  liners  with  very  different 
geometrical  shapes  there  were  not  noticeable  differences  in 
antivibration  effects.  This  shows  that  the  geometric  shape  of 
the  liner  has  a  secondary  influence  on  total  antivibration 
effects . 

(4)  The  antivibration  effects  of  liners  designed  according 
to  the  acoustic  resonance  principle  and  which  are  perforated, 
on  low  frequency  oscillation  is  lacking.  By  increasing  the 
resonant  space  we  can  lower  resonant  space  intrinsic  frequency 
fQ.  Yet  an  excessive  enlargement  of  the  resonant  space  causes 
the  working  condition  speed  in  the  liner's  deep  center  high 
temperature  gas  flow  to  worsen,  for  example  the  type  I  perforated 
acoustic  liner. 

3.  Tests  of  Afterburner's  Shell  Wall  Temperature 

There  are  12  measured  temperature  points  on  the  wall  surface 
of  the  afterburner  shell  I,  II  and  III  sections  to  measure  the 
afterburning  wall  temperature  (see  fig.  1). 

See  figures  11  and  12  for  the  section's  mean  wall  temperature 
changes  along  the  axis  of  the  afterburner  and  the  peripheral 
distribution  of  the  wall  temperatures  in  section  III. 
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Fig.  11  Variation  Curve  of  Mean  Shell  Temperatures  Along  Axis 

Key:  1.  x  (millimeters) 

2.  Type  I  perforated  acoustic  liner 

3.  Type  II  perforated  acoustic  liner 

4.  Type  III  perforated  acoustic  liner 
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12  Peripheral  Distribution  of  Wall  Temperatures  in 
Section  III 

Key:  1.  Type  I  perforated  acoustic  liner 

2.  Type  II  perforated  acoustic  liner 

3.  Type  III  perforated  acoustic  liner 


The  test  results  of  measured  temperatures  show  that: 


(1)  The  enlargement  of  the  secondary  flow  of  the  perforated 
acoustic  liner  and  in  the  shell  crevice  can  cause  the  highest 
mean  wall  temperature  of  the  shell  to  drop.  The  wall  temperature 
of  the  shell ' s  section  III  with  a  type  III  perforated  acoustic 
liner  drops  about  70°C  lower  than  the  shell  with  type  I 
perforated  acoustic  liner.  Its  axial  temperature  rise  gradient 
drops  from  180-200°C/meter  to  110-130°C  meter. 

(2)  A  uniform  secondary  flow  peripheral  distribution  in  the 
crevice  can  cause  the  shell's  peripheral  temperature  difference 
to  noticeably  drop  and  the  peripheral  temperature  difference  of 
the  shell's  section  III  with  a  type  III  perforated  acoustic 
liner  drops  107°C  lower  than  the  shell  with  a  type  I  perforated 
acoustic  liner. 

(3)  The  above  two  factors  cause  the  highest  thermal  point 
temperature  to  noticeably  drop.  The  highest  thermal  point 
temperature  of  section  III  with  a  type  I  perforated  acoustic 
liner  reaches  950°C  and  the  temperature  points  above  900°C 
occupy  one  half  of  the  section.  However,  when  the  highest 
temperature  with  a  type  III  perforated  acoustic  liner  is  only 
860°C,  there  is  a  drop  of  90°c. 

(4)  Analysis  of  Wall  Temperature  Test  Results 

Test  results  show  that  the  geometric  shape  and  measurement 
of  the  perforated  acoustic  liner  has  a  noticeable  influence  on 
the  afterburner's  shell  wall  temperature. 

(1)  The  enlarging  of  the  secondary  flow  of  the  perforated 
acoustic  liner  and  shell  crevice  increased  the  cooling  effect  of 
the  low  temperature  secondary  flow  on  the  snell's  wall  surface. 
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(2)  When  the  liner's  ripple  geometric  shape  was  made  smooth, 
this  caused  the  secondary  flow  distribution  along  the  periphery 
to  be  uniform  which  effectively  weakened  the  high  temperature 
stripes  produced  from  the  backward  stretching  concave  slot 
vortex  at  the  wave  valley  exit  on  the  perforated  acoustic  liner. 

(3)  The  increase  of  the  secondary  flow  in  the  crevice 
decreased  the  amount  of  air  of  the  main  combustion  area  near  the 
perforated  acoustic  liner.  This  caused  the  gas  air  mixture  ratio 
near  the  perforated  acoustic  liner  to  increase,  changed  the 
optimum  chemical  equivalent  ratio  during  combustion,  raised 
combustion  efficiency  in  that  area,  weakened  the  combustion 
elongated  segment  due  to  the  compensation  combustion  of  the 
uncombusted  oil,  and  caused  the  elongated  segment  wall  temperature 
to  rise  slowly. 

IV.  Preliminary  Conclusions 

The  antivibration  effects  of  perforated  acoustic  liners 
mainly  determined  by  the  size  and  perforated  area  ratio  of  the 
resonant  space  and  oscillation  absorptivity  small  holes  in  the 
perforated  acoustic  liner  and  afterburner  shell.  The  perforated 
acoustic  liner's  lowered  effects  on  the  afterburner  shell  are 
not  only  determined  by  the  size  of  the  perforated  acoustic  liner 
and  shell  crevice  but  are  also  determined  to  a  large  extent  by 
the  uniformity  of  the  crevice  along  the  periphery.  This  shows 
that  the  geometric  shape  of  the  perforated  acoustic  liner  is 
only  a  relatively  secondary  factor  for  antivibration  but  it  has 
a  noticeable  influence  on  afterburner  shell  wall  temperature. 
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A  TIME  INVARIANT  STATE  ESTIMATOR  FOR  CONTINUOUS  TIME  SYSTEMS 
by  Dai  Guanzhong 

(Northwestern  Polytechnical  University) 

Abstract 

This  paper  studies  the  design  method  of  the  steady  state 
estimator  for  linear  steady  and  continuous  time  systems  when 
the  system  and  measured  noise  is  a  smooth  white  noise  process. 

To  improve  the  transient  performance  of  the  traditional  Kalman- 
Bucy  filter  definitions  of  two  new  performance  functions  are 
proposed,  thus  giving  two  modified  Kalman-Bucy  filters  which 
can  be  conveniently  applied  in  engineering. 

I.  Preface 

It  is  well  know  that  there  exists  a  basic  contradiction 
between  transient  state  and  steady  state  requirements.  The 
traditional  Kalman-Bucy  filter  has  "narrow  frequency  band" 
problems,  that  is,  there  can  be  an  effectively  filtered  Kalman- 
Bucy  filter  for  the  random  noise  which  is  often  a  relatively 
remarkable  system  of  dynamic  hysteresis  [1-4]. 

When  applied  in  engineering,  the  transient  state  and  steady 
state  performances  of  the  state  estimator  must  be  excellent  and 
one  cannot  be  attended  to  while  losing  sight  of  the  other.  This 
is  especially  the  case  for  state  estimators  operating  for  short 
periods  of  time  where  transient  performance  is  very  important. 

To  improve  the  transient  performance  of  the  traditional 
Kalman-Bucy  filter,  this  paper  proposes  considering  the 
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requirements  of  the  filter's  steady  state  estimated  error  and 
also  defines  two  new  performance  functions  required  for  its 
transient  performance.  In  this  way,  the  paper  gives  the  design 
methods  of  two  modified  Kalman-Bucy  filters  which  can  be 
conveniently  applied  in  engineering. 

II.  System  and  Filter 

The  estimated  linear  steady  and  continuous  time  system  S  is 
described  by  the  following  equations: 

*(  t  )—  Ax(  t  )+Bu(  0  +  r»(  t),  x(0)-x,  (1) 

y  ( t  )=C*(  t )+  V  ( t )  (2) 

In  the  equations,  state  control  u£Rr,  output  y£Rm,  system 

noise  ‘measuring  noise  \*£Rm,-  A,  B,  C  and  T'  are  separately 

the  state  matrix,  control  input  matrix,  output  matrix  and  noise 
imput  matrix  of  (n,n) ,  (n,r),  (m,n)  and  (n,l). 

If  control  £u(t)}  is  the  determined  input,  and  if  £w(t)^ 
and  fv(t)}  are  zero  mean  values,  the  mutually  exclusive  smooth 
white  noise  process  is: 


£w(/)- o,  Ew(i  )»,(t)  =  C?6(  t- t)  (0^0)  (3) 

£»(f)-0,  Ev(t)v'(x)  =  Rb(t  -t)  (R>  0)  (4) 

Eu>(  t  )v'  (*)  =  0  (5) 


Further,  if  initial  state  x  is  an  uncorrelated  random  vector 

o 

with  ^w(t)}  and  fv(t)}  : 
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Ex, w'  (0*0,  Ex,v'  (  t  )  ■■  0 


(«> 

We  know  that  its  mean  value  and  variance  matrix  are: 

Ex,~x„  E(x,-x,)(x,  -x,)'-P,  (/>,>  o)  (7) 

To  guarantee  the  asymptotic  characteristic  requirements  of  the 
filter,  we  also  assume: 

(A,C)  is  a  measureable  pair  (a,  TQ  h)  is  a  controllable 
pair  (8) 

u 

Here,  Q  is  the  square  root  of  the  non-negative  system’s  noise 
variance  matrix  Q,  that  is:  for  Q  >  0,  an  existing  symmetrical 
matrix  M  ^  0  causes  Q=MM*M  ,  noting  M=Q  . 

Because  of  the  corresponding  relation  of  the  primitive  form 
and  model  existing  within  the  estimated  system  and  its  state 
estimator  [5] ,  we  can  assume  that  the  equation  of  the  steady 
state  estimator  (filter)  SE  is: 

i(t)-Fx(i)+L»(l)+Gy(1)'  (9) 

*<1  )-#*(<)  (10) 

In  the  formulas,  the  filter's  state  z£R^,  the  filter's  output 

x£Rn,  which  is  an  estimation  of  state  x;  F,  G,  L  and  H  are 
separately  the  filter's  feedback  matrix,  forward  matrix,  control 
input  matrix  and  output  matrix  of  (p,p),  (p,m) ,  (p,r)  and  (n,p) . 
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Fig.  1  System  and  Filter 

Key:  1.  System  S 
2.  Filter  S£ 

To  sum  up,  the  design  of  filter  S£  is  the  determination  of 

its  dimension  p,  initial  value  z  as  well  as  the  F,G,L  and  H 

o 

matrices. 

It  is  only  necessary  to  accord  with  the  following  state 
estimation  requirements: 

(1)  Unbiasedness,  that  is 

(1) 

a*  £*(O-0  (11) 

Key:  1.  Or 

(2)  Minimum  variance,  that  is 

t,P(  /)—/,£*(  /)*'(  t)  —  min  (12) 

In  the  formulas,  P(t)4  ExCtjx'Ct)  is  the  estimation  error, 

(3)  Transient  performance  requirements,  that  is  the  matching 
requirements  of  the  characteristic  value  (F) (i=l, 2. . . ,p)  for  F. 
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Usually,  there  is  a  contradiction  between  requirements  (2) 
and  (3)  and  it  is  necessary,  based  on  the  different  uses,  to 
lay  stress  on  a  certain  aspect  from  the  different  performance 
functions. 

III.  Unbia  sedne s  s 


From  the  requirements  of  unbiasedness,  we  can  determine 

dimension  p,  initial  value  z  of  St  as  well  as  the  relation 

o 

between  the  F,  G,  L  and  H  matrices. 


Firstly,  when  the  initial  time  t=0,  then  from  the  unbiased¬ 
ness  requirements 

Ex  =Ex  =x 
o  o  o 

A 

Because  the  initial  estimation  X  is  a  nonrandom  quantity, 

o 

therefore 


A  - 
X  — X 

o  o 


(13) 


On  the  other  hand,  the  initial  state  z  of  S_  from  formula  (10) 

o  t, 

should  satisfy: 

„  A  — 

Hz  =  X  =X 
o  o  o 

In  order  for  there  to  be  only  one  solution  of  z  for  an  arbitrary 

x  ,  it  is  then  required  that: 
o 

H€Rnxn(i.e.,  p=n)  and  the  H  matrix  is  non-singular  (14) 

Thus,  formulas  (9)  and  (10)  can  be  changed  into 


t  )-HFH-'Sc(  t)  +  MLu(  i  )  +  HGy(  I  ), 


(15) 


when  comparing  formulas  (1)  and  (15),  in  order  to  causo  the 
influence  of  the  estimated  system  S  and  state  estimator  SE  on 
the  determined  control  u(t)  to  be  the  same,  we  should  take 

HL=B 
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(16) 


Because  of  this,  from  the  t=0  transient  unbiasedness  require¬ 
ments,  the  Sj,  equation  is  then 

-4ri(  t  )-HFH-'H  t)+Bu{t  )+HGy(  t  ),  A,-*,  (17) 

Secondly,  based  on  the  unbiasedness  requirements,  when  t  ^  0, 
we  further  determine  the  relation  -between  the  F, H  and  G  matrices 
in  formula  (17).  We  can  obtain  from  formulas  (1),  (2)  and  (17) 

t  )-Ax(  l)  +  rw(t  t  ) -HG(.Cx(  t)  +  o  (  <)3 

~(A-HFH"-HGC)  x  (  t  )  +  HFH-'Z(  t  )  +  r w(t  )-HGv(  t  ) 

(18) 

From  the  unbiasedness  requirements:  Ex(t)=0,  and  therefore 

dt  ESf(t).EpjtS  (t)J  =0;  further  paying  attention  that  the  mean 

values  of  fw(t)}  and  {v(t)J  are  zero,  we  obtain  from  formula 
(18) 

(A-HFH'-HGC)ExU)- 0  ..  (19) 

On  the  other  hand,  from  formulas  (1)  and  (3)  we  can  obtain 

~^Ex(  t  )-AFx(  t  )+Bu(  t  ),  £*,«%, 

Its  solution  is 

£*(  t )  -«**,+ j  *  *  )dx 
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Generally,  Ex(t)  ^  0,  yet  formula  Cl 9)  is  always  established  for 
t  y  0  and  thus 

A-HFH-1-HGC=0 

For  this  reason,  from  the  unbiasedness  requirements  of  t  ^  0,  we 
can  obtain  the  relation  between  the  F,  H  and  G  matrices  of  S£ 

HFH-'-A-HGC  (20) 

If  formula  (20)  is  substituted  into  formula  (17) ,  we  can 
obtain 

-4fMO-Axit)+Bu(t)+HGiy  (<)-«(.)),  *.-x. 

Forward  matrix  G  and  output  matrix  H  are  joined  into  a  new  forward 
gain  matrix: 

K=HG  (21) 

Then,  we  can  obtain  the  S£  equation 

-^■iOj-^CO+^CO+ifCXO-CiCO),  (22) 

This  type  of  state  estimator  is  unbiased  and  is  called  the 
unbiased  filter.  Fig.  2  shows  the  signal  flow  chart  of  estimated 
system  S  and  unbiased  filter  Sg. 


Fig.  2  System  and  Unbiased  Filter 

Key:  1.  System  S 

2.  Unbiased  filter  S_ 

E 

IV.  Asymptotic  Characteristics 

We  know  from  the  above  discussion  that  for  the  design  of 
filter  SE  the  remaining  requirements  of  (2)  and  (3)  mentioned 
previously  are  used  to  determine  the  forward  gain  matrix  K. 
This  section  first  studies  the  asymptotic  characteristics  of 
filter  SE  and  takes  the  determining  of  the  K  matrix  as 
preparation. 

When  the  unbiased  filter  conditions  of  formulas  (20)  and 
(21)  are  substituted  into  formula  (18) ,  we  can  obtain  an 
equation  for  estimation  error: 

-Jp*(  t  )-(A-KC)x(  t  )  +  Tw(  i  )-Kv(  t  ),  X, =»*„-*,  (23) 

Its  solution  is 

S(  t  )  J '  t  -  t  )  - A'i»(  t  -  t  )Vt  (24) 

When  we  give  attention  to  formulas  (3) -(7),  we  can  obtain  a 
variance  matrix  for  the  estimation  error 

/»<  |  )«#<<*-«®'p,eU-*o/'+  J  '  e '*■**"  (XQT' +KRK,)ei*'KC>,td*  (25) 

Now  we  will  study  the  asymptotic  situation  of  when 

t->a>  in  formula  (25)  ,  if  A-KC  is  a  steady  matrix  then  the 
steady  estimation  error  variance  matrix  is 
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i 


P^limPO)-  f“ +  KRK')et'"u:>''dt  (26) 

[Lemma  1]  If  A-KC  is  a  steady  matrix  and  (a,  p  Q^)  is  a 
controllable  pair,  then  the  P  convergent  of  formula  (25) is 
positive  determine  solved  by  the  following  Lyapunov  equation: 

P(A-KCy  +(A-KC)P  +  (rQrf  +KRK')=  0  (27) 

Proof:  from  formula  (26)  we  can  obtain 


P(A  —KCy  =  e'^'irQr'  +KRK' 

Based  on  partial  integration,  we  obtain 

.  P  ( A-Kcy  =  ror  +  krk*  )*"-kc,',k 

-(A-KC,  J  %' -,-“®*(r<?r/  +  KRK')e  '-“‘"dt 

Taking  notice  that  A-KC  is  a  steady  matrix,  then  from  the  above 
formula  we  know  that  P  satisfies  formula  (27) . 

.  I 

Next,  we  study  the  homogeneous  system  T)  (t)=A-KC)  f)  (t) 
and  take 

V(H  )=  n  Pf) 

Then 


y(r,)*r|' pri+n' Pi=*VC(A-KC)p  +  p(  A-KCy:  n 
-  -  Y  (TQr  +  KRK' )  ’l  «  0 
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Because  A-KC  is  a  steady  matrix,  based  on  the  Lyapunov  steady 
theorem  [6]:  if  the  arbitrary  non  zero  h  Rn,  and  the  V(f))  0 
along  the  solution  locus  >}  (t)=e^A_KC^  '  then  V(T))  is  the 
Lyapunov  function  and  P^O.  In  reality,  if 

v<  t)  )  -  _  *  (ran  +  krk'  )  n  ■  o 


then 


nTOriso,  fOriso 

Because  R  >0,  therefore  K’fjSO  and  also  K=0.  Thus, 

yi 1 P  Q  r  '"1^  =h 1  eAt  P  P  p  q'S  t^_Q 

(for  all  the  t  >  0) 

This  is  a  controllable  relative  contradiction  with  (a ,TQ*).  The 
lemma  proof  is  completed. 

V.  Minimum  Variance 

In  view  of  the  steady  performance  when  the  traditional  steady 
state  Kalman-Bucy  filter  mainly  operates  in  a  t  — *  <*>  steady 
state,  the  performance  function  is  taken  as 

/  +KRK' )ilA'KC>,‘di  (28) 

(Theorem]  As  regards  system  S,  if  (A ,C)  is  a  measurable  pair 
and  when  (a,  V  )  is  a  controllable  pair,  this  causes  the 
performance  function  of  formula  (28)  to  be  minimal  and  the  optimal 
gain  matrix  K  of  filter  SE  is 

K=PC'R_1  (29) 

In  the  formula,  P  is  the  positive  definite  solution  of  the 
following  Riccati  substitution  equation: 
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AP+PA' -PC'R"CP+rQr  -  0 


(30) 


Moreover,  matrix  A-KC  is  a  steady  matrix. 

Proof:  see  reference  [7], 

Formulas  (29)  and  (30)  are  commonly  called  the  substitution 
conditions  of  the  steady  state  Kalman-Bucy  filter  S„.  To  further 

XL 

study  the  transient  performance  of  S_,  we  can  study  the  form  of 

XL 

the  above  mentioned  substitution  conditions  in  the  frequency 
domain.  We  have: 

[Theorem  2]  The  frequency  domain  of  the  steady  state  Kalman- 
Bucy  filter  SE  is 


1/  +R-uiC(jczI-AylKRwi\t 

- 1  +R-1'lc(j<oi-Ay\rQr,)(.-i<oi-A'ylc'R-,'i  oi) 


Here,  R2  is  the  square  root  of  the  positive  definite  noise 

variance  matrix  R,  that  is:  when  R  >  0  and  there  is  a  symmetrical 

/  2  k 

matrix  N  >  0,  R=NN  =N  ,  noting  N=R  . 

Proof:  these  are  also  well  known  results  and  the  proof  can 
be  seen  in  reference  [8]. 

In  particular,  when  estimated  system  S  is  a  single  output 
(m=l)  system,  then  the  frequency  domain  conditions  of  formula 
(31)  in  the  steady  state  Kalman-Bucy  filter  SE  changes  to 


i  nnr' .  jSjamm&e- 
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1 1  +C7  (;W- AY'k\x=  1  +  -€f  (/©/  -  A  )  •*  (ror'  )  ( -  j<oI  -A'Y'c  (32) 


In  the  formula,  S's  output  matrix  c'^R1  .  the  measured  noise's 

2 

variance  R=r  (scalar  quantity),  and  S's  optimal  gain  matrix 
lc£RnXl . 


-1, 


We  know  form  fig.  2  that  G (s) =c 1 (sI-A) _ ^is  the  open  loop 
transfer  function  of  filter  S£  and  therefore  the  .frequency 
domain  conditions  of  formula  (32)  changed  to 


i  +  g  (ym)!1*  i  +-jrc,(j«>i-AY'(rQn)(-j<t>i-A,y'c 


(33) 


that  is 


J1+C0®)|>1  (34) 

This  signifies  that  the  distance  between  the  open  loop  Nyquist 
curve  G(jtf>)  of  the  steady  state  Kalman-Bucy  filter  and 
(-1 , j 0)  point  on  the  complex  plane  is  not  smaller  than  1.  In 
other  words,  G(j©)  cannot  be  advanced  into  the  unit  circle 
taking  the  ( — 1 , j 0 )  as  the  center  of  the  circle  (see  fig.  3). 
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Fig.  3  Restricted  Zone  of  G(j<5>) 

Key:  1.  Restricted  zone 

It  is  easily  known  that  the  phase  capacity  of  filter  S_  is  not 
o  k 

less  than  60  and  the  gain  capacity  is  unlimited  [8] . 

Although  the  steady  state  Kalman-Bucy  filter  S£  has  relatively 
good  asymptotic  steady  characteristics  yet  its  transient 
performance  is  not  good  enough.  In  reality,  S£  is  only  a  "narrow 
frequency  band"  low  pass  filter. 

VI.  One  Method  For  Improving  Transient  Performance 
We  define  the  new  performance  function  as 

/.  -  UP.  -  U  j  ~  eI"Ce  '-kc"(rQr'  +  KRK'  )e^KC> "]  di  (35) 

In  the  formula,  o>0,  and  the  dimension  is  sec-1.  Formula  (35) 
can  be  changed  to 
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/.  =  /,P.  =  ;,  j  ~  eiA*m,'KC>'( TQr'  +KRK'  )e'A*’,'KC,'id{ 


(36) 


In  formula  (36) ,  the  state  matrix  is  already  changed  to  A+al 
and  the  remainder  are  the  same  as  performance  function  (28) . 

[Lemma  2]  (A+a/,,C)  is  the  measureable  pair  but  only  (A,C)  acts 

as  the  measureable  pair,  and  (A+ai,.  V  Q^)  is  the  controllable 
pair  but  only  (A,Tq2)  acts  as  the  controllable  pair. 

Proof:  take  note  that  the  measureable  matrix  is 


c 

1 

r  c  i 

CA 

C(AAal) 

.  CAT1 

j 

lC  ( .1  +«/)“**. 

In  the  formulas 


C  (  A  +  a  / )  —  CA  +  aC 
C  (  A  +  a/  y-  =  CA*  -r  2aC.  J  ru:C 


Because  after  the  constant  of  a  certain  line  multiplies  to 
another  line  the  rank  of  the  matrix  is  unchanged,  therefore 

rank  R  =rank  u 

O 

Proof  for  the  front  half  of  the  lemma  is  completed. 

For  the  back  half  of  the  lemma  we  <-c.'  use  the  duality 
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principle  proof  of  controllability  and  measureability . 

Actually,  from^  the  duality  principle,  (a,  V  Q*5)  is  a  controllable 

pair  and  (a/  , o'5  P  )  is  a  measureable  pair;  from  the  front  half  of 

this  lemma,  [  (A  +<*/)',  Q  ^  T  ’  ]  is  a  measureable  pair  and  (A+d,  fQ^) 
is  a  controllable  pair.  The  lemma  proof  is  completed. 

After  attaining  lemma  2,  the  entire  results  in  sections  IV 
and  V  can  be  copied. 

Firstly,  based  on  lemma  1,  if  A4a/-KC  is  a  st^eady  matrix 
then  the  P  convergence  in  formula  (36)  is  positive  definite 
solved  by  the  following  Lyapunov  equation: 

P.(.A+*I~KCy  +  (A+aI-KC)  P  +  (rQn  +KRK')=  0  (37) 

Secondly,  based  on  theorem  1,  we  can  cause  formula  (36) 
which  is  the  smallest  optimal  gain  matrix  j(m  to  be 

K.~P.C'R “  (38) 

In  the  formula,  pw  is  the  positive  definite  solution  of  the 
following  Riccati  substitution  equation: 

(A+«i)P.+p.(A+«iy -p.c'  R-lcp.+rQr'  =  o  (39) 

Moreover,  matrix  a + a/- K.C  is  a  steady  matrix. 

To  sum  up,  we  then  have  the  following  theorem: 

[Theorem  3]  For  system  S,  if  (A,C)  is  a  measureable  pair  and 
(a,  r  q*5;  is  a  controllable  pair,  then  the  optimal  gain  matrix 
K.  of  the  modified  steady  state  Kalman-Bucy  filter  S.  satisfies 
formulas  (38)  and  (39) . 


We  will  now  analyze  the  transient  performance  of  the  modified 
Kalman-Bucy  filter  S. 

[Theorem  4]  The  distance  of  the  characteristic  value 
distribution  of  the  modified  Kalman-Bucy  filter  S.  in  the  left 
half  of  the  complex  plane  and  of  the  imaginary  distance  axis 
is  larger  than  a. 

Proof:  because  matrix a-¥<U—KjC  is  a  steady  matrix,  its 
characteristic  value  is  distributed  in  the  left  side  of  the 
complex  plane.  On  the  other  hand,  the  corresponding  characteristic 
values  of  matrix  A—K.C  characteristic  value  relative  matrix 
A  +aI—K.C  shift  a  distance  of  a  to  the  left.  Thus,  this  theorem 
is  completed. 

we  saw  that  the  attenuation  of  filter  S.  is  greater  than  a, 
in  other  words,  the  time  constant  of  the  major  extreme  point 
of  s.  is  smaller  than  l/<*. 

We  will  now  further  compare  the  transient  performances  of 
Kalman-Bucy  filter  S£  and  modified  Kalman-Bucy  filter  S’.  For 
this  reason,  it  is  the  same  as  theorem  2  and  thus  we  have: 

[Theorem  5]  The  frequency  domain  conditions  of  modified 
steady  state  Kalman-Bucy  filter  S.  is 
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T 


1/  +  R'l'*C(i<oI -  A  ) "AT.# J/Y-  /  +  Rm,/tC(](itI 

-  a  y\ror  +  zap.)  ( -  ;w  -  a'  ylc'  r-1'*  (40) 


Proof :  We  use  the  Riccati  substitution  equation  to  rewrite 
formula  (39)  as 

AP.  +  P.A' -P.C'R-‘CP.  + (rQr' +  2aP.)  =>  o  '  (41) 

Further  imitating  theorem  2  we  can  immediately  prove  that  the 
frequency  domain  conditions  of  S  are  formula  (40) . 

St 

When  estimated  system  S  is  a  single  output  system,  the  -S', 
frequency  domain  conditions  of  (40)  change  to 

1 1  +✓</«>/- .4 1  +-p-c'(ja>I- A)  '(rQn  +2aP.)(-j<oI-A'y'c  <42) 


If  G.(  *  )~c/ (si  —  Ay'k.  is  the  open  loop  transfer  function  of 

filter  S.  and  we  also  take  into  account  formula  (33) ,  then 
formula  (42)  can  be  rewritten  as 

|1  +<?.0«D)|*-|l  +  G(/<a)|t+-p  cf  A  )’]  P.(  —  j<*>I  —  A'  )''c 

Because  a>o  and  P.  are  positive  definite,  therefore  aside  from 
the  &>  values,  there  is 

•  |l  +C  (;*»)! 
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(48) 


This  signifies  that  the  open  loop  Nyquist  curve  <7.(/<d)  of  St 


is  further  from  the  ( —1 , j 0 )  point  on  the  complex  plane  than 
the  open  loop  Nyquist  curve  G0<°)  of  Sg.  It  is  commonly  known 
that  generally  speaking  the  open  loop  Nyquist  curve  of  the  unit 
feedback  system  is  further  from  the  (-l,jO)  point  and  its  close 
loop  major  extreme  point  is  further  from  the  imaginary  axis. 
Therefore,  the  modified  steady  state  Kalman-Bucy  filter 
raised  transient  performance. 

VII.  Second  Method  For  Improving  Transient  Performance 

Formula  (25)  gave  the  changes  of  the  estimation  error 
variance  matrix.  The  first  item  is  the  transfer  of  the  estimation 
error  caused  by  the  initial  estimation  error  variance  matrix  P  , 
and  the  second  item  is  the  accumulation  of  the  estimation  error 
variance  matrix  produced  by  system  noise  £w(t)^  and  measured 
noise  fv(t)^  .  The  transient  performance  of  the  filter  can  be 
determined  by  the  attenuation  speed  of  the  initial  estimation 
error  Variance  matrix.  For  this  reason,  we  defined  another  new 
performance  function 

(ran  +KRKf')*'‘-™''dt+v  J*  •'--«>"<// } 

(44) 

In  the  formula,  /ff  ^  0.  Here,  the  first  item  is  the  accumulated 
steady  state  estimation  errcr  variance  matrix  produced  by  the 
system  noise  and  measured  noise.  It  delineated  the  steady  state 
precision  of  the  filter.  The  second  item  represents  the  total 
accumulation  effect  produced  by  the  initial  estimation  error 
variance  matrix.  It  delineates  the  transient  performance  of  the 
filter  in  a  fixed  range.  $  is  then  the  weight  coefficient 
between  the  steady  state  precision  requirements  and  transient 
performance  requirements.  Therefore,  can  be  called  the 

comprehensive  steady  state  estimation  error  variance  matrix. 
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Formula  (44)  can  be  rewritten  as 


*<*mKO,(8Pt  +  rQr'  +KRK')e1*-**'’  di  (45) 

In  the  formula,  the  non-negative  Q  PQ+  VqT  '  can  be  decomposed 

into:  /ffPQ=  f1  Q  T  '  ='MN  and  M  is  the  non-negative  symmetrical 
matrix. 

Therefore ,  we  can  have : 

[Lemma  3]  If  (a,  V  Q^)  is  a  controllable  pair  then  (A,M) 
must  be  a  controllable  pair. 

Proof:  If  (A,  V  Q^)  is  a  controllable  pair,  then  it  must 
exist  in  T  <  co  causing 

j  *  e*TQut  (rQ11)'  e*"  di  -  J  ^  e*TQr'  e*"  di>  0 


Further,  because  /?>0,  the  Pq  matrix  is  non-negative  and 
therefore 


J  MM'  e di  -  J  *  (&P,  +  rQr )  H"  di 

-  e"  rQr'eJ,'d1+V  e''Pydi>* 

This  proves  that  (A,M)  is  a  controllable  pair. 

After  obtaining  lemma  3  the  entire  results  of  sections  IV 
and  V  can  be  copied. 

Firstly,  based  on  lemma  1,  if  A-KC  is  a  steady  matrix  then 
the  P^>  congruence  in  formula  (45)  can  be  positive  definite 
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solved  by  the  following  Lyapunov  equation: 


P„(A-KCy  +  (A-KC)Pll  +  (&Pt  +  rQr'  +KKK')  =  0  (46) 

Secondly,  we  also  have: 

[Theorem  6]  For  system  S,  if  (A,C)  is  a  measureable  pair, 
when  (k,  is  a  controllable  pair  then  the  performance 

function  of  formula  (45)  is  minimal  and  the  optimal  gain  matrix 
of  the  modified  steady  state  Kalman-Bucy  filter  sa  is 

Kt=P»C'R-'  '  (47) 

In  the  formula,  Pg  is  the  positive  definite  solution  of  the 
following  Riccati  substitution  equation: 

AP'  +  PtA'  -P,C' R-'CP,  +  (&Pt+ TQr' )~  0  (48) 

Moreover,  matrix  A-I^C  is  a  steady  matrix.  The  frequency  domain 
condition  of  is 


|  /  +  /?*,'*C(;W-  AY  KtR'  l\'-  I  +/?'1/1C(/a>/ 

-  A)‘(0Pt+ rQr  )(-/©/- a  y'C'R-''*  (49) 


Proof:  proof  obtained  from  theorems  1  and  2. 

To  compare  the  traditional  Kalman-Bucy  filter  SE  with  /Q  =0 
and  the  transient  performance  of  the  modified  Kalman-Bucy  filter 
S g  with  /&  y  0 ,  we  studied  the  situation  wherein  estimated  S  was 
a  single  output  system.  At  this  time,  the  filter  frequency 
domain  conditions  of  formula  (49)  changed  to 


li  +tf'ow-yO-%i‘=.  i  +-p-c'u<oi-Ayt(ppt+rQr,)(-j<oi-A'yic(s9) 

If  G^  (s}=c  (sI-A)  is  the  open  loop  transfer  function  of 

filter  S £  and  we  consider  formula  (33) ,  then  formula  (50) 
can  be  rewritten  as 

|1  +<7,(;o>)|1=|  l  +<?0'oi>)|*+-p— s' A)~iPt(  —  jci>I—A>y,c 
Because  (3  y 0  and  Pq  are  non-negative,  therefore 
jl  +<J»(;flo)|>|  1  +G(jffl)J  (51) 

.  r  > 

This  signifies  that  the  open  loop  Nyquist  curve  (j#)  of  Sg 
is  further  from  the  ( —1 , j 0 )  point  on  the  complex  plane  than  the 
open  loop  Nyquist  curve  G(j<3>  )  of  Sg.  Therefore,  the  modified 
steady  state  Kalman-Bucy  filter  raised  transient  performance, 

VIII.  Conclusion 

Based  on  the  different  requirements  for  application  in 
engineering,  the  designs  for  the  steady  state  estimator  for 
continuous  time  systems  can  be  carried  out  in  three  cases: 

(1)  Transient  design  in  view  of  the  transient  performance. 
This  is  the  design  problem  of  a  state  observer; 

(2)  Steady  state  design  in  view  of  steady  state  precision. 
This  is  the  design  problem  of  the  traditional  Kalman-Bucy 
filter? 
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PW©r«:"~-L.i(*v 


►  V  -  :  •a*':- 


(3)  When  considering  the  steady  state  and  transient  state 


synthetic  design  with  steady  state  precision  and  transient 
performance,  this  paper  proposed  design  methods  for  two  types 
of  modified  Kalman-Bucy  filters. 
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SELF-EXCITED  RESONANCE  OF  AEROSTATIC  BEARINGS  IN  A  VACUUM 
ENVIRONMENT 

by  Fu  Xianluo 

(Institute  of  Mechanics,  Academic  Sinica) 

Abstract 

The  cause  of  instabilities  in  an  aerostatic  bearing  system 
operated  in  a  vacuum  environment  is  studied  in  this  paper.  In 
order  to  eliminate  this  type  of  self-excited  resonance,  some 
suitable  measures  are  suggested  which  can  be  used  to  maintain 
the  ambient  pressure  of  the  bearing  in  a  range  of  about  ^  0.2 
atm. .  This  will  improve  the  precision  of  the  instruments  using 
aerostatic  bearings  as  well  as  increase  the  percentage  yield 
of  qualified  bearings. 

This  work  was  done  under  the  guidance  of  comrade  Lin  Tongji. 

I.  Suggestion  of  Problem 

The  application  of  aerostatic  bearings  has  become  increasingly 
widespread  in  many  areas  of  industry,  scientific  research  and 
national  defense.  Because  their  frictional  force  is  small,  they 
possess  the  special  features  of  small  power  consumption  and  cool 
state  operation;  they  can  operate  at  an  extremely  high  rotating 
speed,  their  axis  rotating  precision  is  very  high  and  abrasion 
wear  extremely  small  and  therefore  has  a  long  operating  life; 
they  do  not  very  much  require  or  basically  do  not  require 
scheduled  maintenance  but  can  operate  in  extremely  high  or 
extremely  low  temperature  environments;  their  noise  and  vibration 
are  small.  Therefore  the  application  of  aerostatic  bearings  is 
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especially  widespread  in  precision  machine  tools  and  precision 
instruments . 

However,  it  is  not  that  easy  for  the  bearings  to  realize  the 
above  mentioned  special  features  under  special  application 
situations.  For  example,  when  the  ambient  pressure  of  the 
bearing  is  very  low,  in  a  high  vacuum  environment,  aerostatic 
bearings  can  often  produce  self-excited  resonance  instability 
problems.  The  higher  the  environmental  vacuum  the  easier  it  is 
for  the  bearing  to  produce  self-excited  resonance.  To  avoid  this 
type  of  resonance,  bearings  measurements  must  be' designed 
accurately  and  their  machining  precision  requirements  are  also 
very  high.  Therefore,  when  aerostatic  bearings  are  used  in  a 
high  vacuum  environment,  the  slightest  change  in  bearing 
measurements  can  produce  self -excited  resonance  and  thus  destroy 
the  normal  operation  of  the  precision  instruments.  Furthermore, 
the  self-excited  resonance  of  aerostatic  bearings  in  this  high 
environment  vacuum  is  also  very  sensitive  to  other  factors. 

Errors  during  machining  and  measuring,  changes  of  the 
environmental  temperature,  the  intrinsic  frequency  of  bearings, 
mounting  methods  and  air  feed  pressure  can  possibly  bring  about 
self-excited  resonance  of  bearings  and  cause  the  normal  operation 
of  the  aerostatic  bearings  to  break  down. 

II.  Qualitative  Analysis 

During  the  design  of  any  aerostatic  bearing,  the  guarantee 
that  it  can  rotate  steadily  under  certain  conditions  is  an 
important  factor  to  consider.  Naturally,  if  the  bearing  system 
is  unstable  then  before  eliminating  the  instability  we  consider 
that  the  other  aspects  of  the  design  are  worthless.  If  we  do  not 
consider  the  problem  of  stability  during  design  or  the 
considerations  are  incomplete  then  the  vast  majority  of 
aerostatic  bearings  can  manifest  self-excited  resonance. 


144 


It  has  been  practically  proven  that  when  the  aerostatic  bearings 
of  a  precision  instrument  is  in  a  high  vacuum  environment  it  can 
manifest  self-excited  resonance.  During  resonance  a  sound  can  be 
heard  and  the  resonance  can  also  be  felt  with  the  hand,  When 
serious,  the  bearing  makes  a  shrill  sound  which  shifts  the  entire 
base. 

The  problems  of  the  instabilities  in  aerostatic  bearings  has 
attracted  the  attention  of  many  researchers  who  have  carried  out 
theoretical  analyses  of  bearings  with  simple  geometric  shapes 
and  have  adopted  stabilization  measures  [1-3].  However,  among 
many  of  the  actual  problems,  because  the  geometric  shapes  of  the 
bearings  were  complex  (for  example,  some  aerostatic  bearings 
including  journal  bearings  and  the  Yates  bearings  of  thrust 
bearings) ,  because  it  was  necessary  to  use  several  of  these  types 
of  bearings  in  the  instrument  and  because  the  mounting  positions 
were  not  matched,  the  temperature  changes  and  loads  endured  by 
the  bearings  was  very  large.  This  caused  the  calculations  of  the 
aerostatic  bearings’  instabilities  to  be  extremely  difficult. 

Even  if  the  geometric  shape  were  simplified  and  ideal  operation 
conditions  were  used  so  that  the  set  of  equations  for  solving 
the  fluid  kinetics  becomes  possible,  we  must  still  use  an 
electronic  computer  to  carry  out  long  computations.  Done  in  this 
way,  the  cost  is  very  high  and  the  obtained  results  can  be  quite 
different  from  the  test  results. 

Therefore,  we  do  not  follow  the  previously  used  method  but 
simplify  the  complex  geometric  shape  and  appropriately  idealize 
the  severe  operation  conditions  so  that  it  is  possible  to  solve 
the  set  of  kinetic  equations.  We  attempted  to  carry  out 
qualitative  analysis  of  actual  phenomena.  Thus,  we  found  the 
source  of  the  resonance  and  adopted  suitable  measures  to  change 
the  position  of  the  resonance  source,  and  weaken  or  even  eliminate 
the  resonance  source. 
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As  regards  the  aerostatic  bearing  shown  in  fig.  1,  when 
the  ambient  pressure  is  about  atmospheric  pressure,  operation 
is  normal.  When  the  environmental  vacuum  is  very  high,  self- 
excited  resonance  can  appear. 


Fig.  1  Scheme  of  the  Aerostatic  Bearing 

Based  on  tests,  it  is  only  necessary  to  continually  raise 
the  environmental  vacuum  for  this  type  of  aerostatic  bearing  to 
produce  self-excited  resonance.  Given  this  point,  we  can  infer 
that  the  resonance  source  that  brought  about  the  self -excited 
resonance  is  possibly  a  supersonic  flow  appearing  in  a  certain 
area  of  the  bearing.  Naturally,  a  lowered  resonance  pressure  P 

s 

can  cause  the  vacuity  of  the  self -excited  resonance  environment 
to  rise.  Yet,  lowered  resonance  pressure  directly  influenced 
the  decrease  of  the  lubricated  gas  flow  and  as  a  result  the 
bearing  capacity  of  the  bearing  dropped.  Although  this  raised 
the  resonant  vacuity  of  the  self-excited  resonance  we  were  not 
able  to  estimate  the  maximum  load  of  the  instrument  during 
design . 

We  know  from  the  theories  of  aerodynamics  that  to  produce  a 
supersonic  air  flow  for  a  one  dimensional  flow,  as  regards 
geometric  shape,  it  is  only  necessary  to  have  a  first  contraction 
then  expansion  passage.  We  can  see  from  fig.  1  that  only  if  there 
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is  an  annular  intake  narrow  crevice  AB  and  throttling  step  CD 
area  on  the  'hrust  surface  can  there  be  air  flow  stoppage. 

If  the  resonance  source  is  formed  from  air  flow  stoppage 
caused  by  the  annular  intake  narrow  crevice  AB  area,  then 
pressures  P0  in  the  B  point  area  must  satisfy 

PB/Ps  ^  0-528  ('>'=1*4) 

Because  Pg  is  positive  definite,  in  order  for  narrow  crevice  AB 

not  to  become  the  source  of  resonance,  it  is  necessary  to  raise 

the  pressure  P  at  the  B  point,  rendering 
B 

P_/P  >0.528 

B  S 

When  the  width  of  narrow  crevice  AB  is  Spu  this  produces 

self-excited  resonance  and  the  width  of  the  narrow  crevice 

increases  to  (5  +1)/*,  •  When  the  air  flow  passes  narrow  crevice 

AB  the  loss  will  decrease  along  the  way.  Thus,  P  is  larger  than 

B 

the  original  value  and  it  possibly  does  not  form  the  supersonic 
resonance  source.  Nevertheless,  tests  have  proven  that  by 
enlarging  the  width  of  narrow  crevice  AB,  self-excited  resonance 
appears  at  even  lower  vacuity  than  the  original  vacuity  and 
therefore  narrow  crevice  AB  is  not  the  source  of  resonance. 
Because  the  possible  source  of  resonance  is  only  in  two  areas 
the  source  of  resonance  can  only  be  throttling  step  CD. 

We  can  also  analyze  that  the  source  of  resonance  is  the 
throttling  step  CD  from  two  other  aspects.  Firstly,  the  enlarging 
of  the  width  of  narrow  crevice  AB  caused  resonating  vacuity  to 
drop.  We  can  see  that  the  enlargement  of  the  intake  narrow 
crevice  decreased  the  loss  along  the  way  and  thus  raised  P  as 

B 

well  as  Pfi.  This  can  possibly  cause  pD/pc  ^  0.528.  In  this 

way  a  supersonic  air  flow  is  formed  beginning  from  the  back  edge 
of  throttling  step.  Secondly,  based  on  tests,  when  we  eliminated 
the  thrust  end  cover  of  the  bearing,  the  resonating  environment's 
vacuity  noticeably  rose.  It  can  also  be  concluded  that  throttling 


Fig.  2  Schematic  of  Pressure  Distribution 

We  will  use  the  pressure  distribution  schematic  shown  in 

fig.  2  to  further  explain.  If  the  air  flow  in  the  bearing  is 

an  axial  symmetrical  one  dimensional  flow,  the  pressure  is  a 

linear  distribution,  x  is  the  flow  line  coordinate,  the 

resonance  pressure  in  the  A  point  area  is  P  ,  and  the  ambient 

s 

pressure  is  P  ,  the  pressures  of  points,  B,C  and  D  are  P  ,  P 
a  d  c 

and  PD> 

If  pressure  distribution  curve  a  corresponds  to  the  width  of 

the  intake  narrow  crevice  which  is  Sjm  .  under  these  conditions, 

when  ambient  pressure  is  P  self-excited  resonance  is  produced 

a 

and  naturally  Pa/pc  ^  0.528.  If  the  width  of  intake  narrow 

crevice  AB  is  decreased,  for  example  decreased  to  (  6  -l  )/*.  then 

P  drops  and  possibly  causes  P  /P  >  0.528  and  the  air  flow  in 
c  a  c 

front  and  behind  the  throttling  step  to  be  a  subsonic  flow.  This 

corresponds  to  the  situation  of  pressure  distribution  curve  B. 

If  environmental  vacuity  is  continually  raised  wherein  the  value 

of  P  drops,  this  causes  the  value  of  P  /P  to  drop  resulting  in 
a  a  c 


a  supersonic  flow.  If,  at  this  time,  the  width  of  the  narrow 
crevice  is  continually  decreased,  for  example  it  is  ($  -2 )///, 
the  value  of  P  /P  corresponding  to  curve  Y  also  increases,  the 

9  C 

air  flow  becomes  a  subsonic  flow  and  the  resonance  source  is 
eliminated . 

We  can  see  that  if  the  environmental  vacuity  is  raised 
(P  is  lowered) ,  we  can  decrease  the  width  of  the  intake  narrow 
crevice  to  eliminate  the  source  of  resonance.  However,  if  the 
narrow  crevice  continually  narrows  the  rate  of  flow  within  the 
bearing  will  become  smaller  and  smaller  and  the  bearing  capacity 
of  the  bearings  will  also  become  increasingly  smaller  and  not 
accord  with  design  requirements. 

III.  Orifice  Throttling 

We  can  see  from  the  above  qualitative  analysis  that  if 
ambient  pressure  continually  drops  and  at  the  same  time  a  fixed 
bearing  capacity  is  guaranteed,  then  the  self-excited  resonance 
of  the  bearing  is  unavoidable. 

To  eliminate  the  resonance  source  on  an  aerostatic  bearing 
and  guarantee  that  the  bearing  itself  will  not  produce  self¬ 
exciting  resonance  under  any  ambient  pressure  but  will  be  able 
to  operate  normally,  we  can  use  a  sealing  cover  to  seal  up  the 
entire  instrument  using  aerostatic  bearings.  At  the  same  time,  we 
process  a  critical  nozzle  (open  a  specially  designated  orifice) 
on  the  shell  of  this  sealing  cover.  After  causing  the  ambient 
pressure  to  be  lower  than  a  certain  value,  even  if  the  ambient 
pressure  continually  drops,  the  pressure  in  the  sealing  cover 
can  be  maintained  lower  than  (or  close  to)  a  certain  design 
value  of  atmospheric  pressure.  Furthermore,  for  this  type  of 
critical  nozzle,  when  the  ambient  pressure  is  an  atmospheric 
pressure,  it  is  also  necessary  to  cause  the  pressure  in  the 
sealing  cover  to  be  greater  than  (or  close  to)  a  certain  design 


value  of  an  atmospheric  pressure.  Done  in  this  way,  the 
resonance  source  in  the  bearing  is  shifted  to  the  shell  of  the 
sealing  cover. 

To  find  the  size  of  the  critical  nozzle's  aperture  we 
analyzed  this  problem  according  to  the  one  dimensional  flow 
theory . 


The  equation  of  conservation  of  momentum  for  ideal  fluid  one 
dimensional  constant  flows  is  [4] 


du  _ 1 _ dp_ 

dx  "  P  dx 


(1) 


^uA=a  constant  (2) 

Based  on  the  definition  of  the  speed  of  sound,  we  have 


we  can  change  equation  (1)  into 


- (3) 

The  two  sides  of  equation  (2)  are  logarithmic  integrals  and 
when  the  results  are  joined  with  equation  (3)  and  we  eliminate 
d  f>  /  P  ,  we  obtain 

du  _  1  dA  /,) 

u  M’-l  A 


In  the  formula 


M— i“  (5) 

is  the  Mach  number.  We  took  the  logarithmic  integral  for  equation 
(5)  and  for  the  Bernoulli  equation 

2  2 

^2  +  =  a  constant 
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integral,  we  can  obtain 

TT-(1+— 

t 

Comparing  the  integrals  of  equations  (4)  and  (.6)  ,  we  have  [4] 

y  *  j 


(7) 


In  the  formula,  A*  is  the  critical  section  of  the  critical  nozzle. 

Because  the  volume  of  the  sealing  cover  is  large  we  can  take  the 

flow  velocity  of  the  gas  in  the  cover  and  the  Mach  number  as 

being  equal  to  zero.  We  used  p  ,  P  and  T  to  indicate  the 

o  o  o 

thermodynamic  parameters  of  the  gas  in  the  sealing  cover.  If  the 
exit  section  area  of  the  orifice  is  A  and  counterpressure 
(ambient  pressure)  is  p',  we  can  deduce  that  the  per  second  mass 
flow  m  of  the  gas  is 


m  =* 


Y  »  1 
2  (  Y  -  l  ) 


v'y P,P,  AW 


Maximum  flow  m*  is 


(  8  ) 


Y  <•  1 


Syp.pt  a 


(9) 


If  the  M  number  in  function  .0  is  indicated  as  the  function  of 

pVp_/  the  [4] 
o 
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1 

I 

1 

1 

) 

_ _ _ _ _  i 

M- <>•> 


If  m=m*,  that  is,  when  the  outlet  velocity  of  the  gas  in  the 

orifice  on  the  sealing  cover's  shell  reaches  the  'speed  of 

sound,  the  p'/p  =0.528  (taking  Y=1.4v  If  pressure  p  in  the 
o  )  .  o 

cover  must  be  less  than  0.84  atmospheric  pressure,  then 

counterpressure  p'  should  be  less  than  or  equal  to  0.44 

atmospheric  pressure.  In  other  words,  when  environmental  vacuity 

is  equal  to  or  higher  than  this  value,  the  pressure  in  the  cover 

is  maintained  at  0.84  atmospheric  pressure. 

Assuming  that  the  gas  mass  flow  flowing  into  the  cover  is  Q 
grams/second  and  the  outlet  speed  of  the  orifice  on  its  shell 
reaches  the  local  speed  of  sound,  then  the  radius  of  the  orifice 
is 


R 


/ 


Qt 


Y  +  1 


✓  Y  P% 


(U) 


If  the  flow  flowing  into  the  sealing  cover  is  invariant  and 
the  size  of  the  orifice  is  also  invariant  but  the  counterpressure 
rises  to  an  atmospheric  pressure,  then  the  pressure  of  the  gas 
in  the  sealing  cover  should  satisfy  the  following  formula: 
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Under  our  conditions,  when  taking  V=1.4,  the  above  formula  can 
be  changed  to: 


(12) 


In  order  to  differentiate,  we  used  pQ^  and  in  the  above 

formula  to  indicate  that  the  orifice's  outlet  speed  is  the  local 
speed  of  sound,  and  the  counterpressure  is  the  pressure  in  the 
sealing  cover  when  there  is  atmospheric  pressure. 


t  ♦  1 


To  make  the  selection  by  designers  convenient,  table  1  lists 
the  highest  counterpressures  of  supersonic  flows  brought  about 
under  specific  flow<?=1.5  grams/second  conditions  as  well  as  the 
corresponding  pressures  pQ^  and  pq2  orifice  critical  sectional 
area  and  its  radius  R. 

It  should  be  pointed  out  that  when  the  ambient  pressure  is 
an  atmospheric  pressure,  due  to  the  influence  of  the  boundary 
layer  formed  in  the  throttling  orifice,  the  effective  area  of 
the  small  hole  shrinks  and  therefore  pressure  pQ2  in  the  cover 
is  slightly  larger  than  the  values  listed  in  table  1. 
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Table  1  Parameters  Associated  With  Throttling  Through  An 
Orifice 

Key:  1.  (Atmospheric  pressure) 

2.  (Atmospheric  pressure) 

3.  (Atmospheric  pressure) 

2 

4.  (Centimeters  ) 

5.  (Centimerers) 

The  instruments  with  aerostatic  bearings  are  placed  in  the 
cover  and  additional  throttling  orifices  are  opened  on  the  cover's 
shell.  This  raises  the  pressure  in  the  cover  and  eliminates  the 
resonance  of  the  bearing  itself.  However,  resonance  caused  by 
supersonic  flow  appears  in  the  outlet  area  of  the  throttling 
orifices  on  the  sealing  cover's  shell.  The  smaller  the  orifices 
the  higher  the  resonance  frequency  and  the  larger  the  orifice's 
outlet  speed  the  greater  the  resonating  energy.  Actually,  the 
method  of  throttling  through  orifices  is  to  shift  the  self-excited 
resonance  of  the  bearing  itself  to  the  sealing  cover's  shell. 

Shifting  the  resonance  of  the  bearing  in  the  sealing  cover 
to  the  cover's  shell  is  noticeably  beneficial  because  the 
resonance  of  the  bearing  causes  precision  instruments  to  work 
improperly.  The  resonance  of  the  sealing  cover's  shell  only 
influences  the  precision  of  instruments. 

Based  on  experiences,  the  influence  of  the  resonance  caused 
by  the  orifice  throttling  in  the  sealed  cover's  shell  on 
instrument  precision  can  be  overlooked  but  we  should  use  this 
test  to  measure  the  level  of  this  type  of  influence. 

If  after  measurements  this  type  of  influence  is  found  to  be 
formidable,  then  there  are  also  means  to  eliminate  the  influence 
of  this  type  of  resonance  on  the  instruments.  For  example,  we 
can  use  a  special  tube  which  can  easily  transfer  resonance 
(its  inner  diameter  must  be  much  larger  than  the  orifice's 
diameter).  The  sealing  cover's  shell  is  joined  to  the  throttling 


155 


orifice  and  damping  equipment  is  used  to  fix  the  throttling 
orifice  on  the  shell. 

We  can  also  directly  open  the  throttling  orifice  on  the 
shell  and  mount  a  muffler  in  the  downstream  of  the  orifice. 

The  muffler  can  be  designed  according  to  information  given  in 
related  references  [5] .  However,  we  should  pay  attention  that 
the  muffler  is  designed  carefully  otherwise  there  is  the 
possibility  that  it  will  be  unable  to  muffle  sound  but  rather 
will  cause  resonance  and  intensify  the  resonance  of  the  sealing 
cover's  shell. 

Fine  porous  throttling  can  also  decrease  resonance,  that  is, 
a  single  throttling  orifice  is  made  into  many  fine  orifices. 

When  the  air  flow  flows  in  the  fine  orifices  the  viscous  formed 
boundary  layer  decreases the  effective  area  of  the  fine  orifices. 
Therefore,  the  area  of  the  fine  orifices  should  be  larger  than 
the  area  of  a  single  orifice.  The  finer  the  fine  orifices  the 
larger  their  area.  However,  they  cannot  exceed  three  times  the 
area  of  an  orifice  when  there  is  single  orifice  throttling  (this 
value  will  be  explained  in  the  next  section  "seepage  flow 
throttling") .  Concrete  details  can  be  determined  from  tests. 

Seepage- flow  throttling  discussed  in  the  next  section  is  also 
a  method  which  can  be  used  for  selection. 

IV.  Seepage  Flow  Throttling 

As  regards  fine  porous  throttling,  if  the  orifice's  diameter 
is  selected  on  an  increasingly  finer  basis  then  the  viscosity  of 
the  gas  becomes  the  dominant  effect.  When  the  orifice's  diameter 
is  very  fine,  which  even  machining  cannot  bring  about,  we  can  use 
a  porous  metal.  This  is  seepage  flow  throttling.  Because  the  flow 
speed  of  seepage  flow  is  very  small,  seepage  flow  throttling 
cannot  manifest  supersonic  resonance. 
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Based  on  the  Darcy  law  [6,7]  of  the  seepage  flow  theory,  we 
can  derive 


m 


(14) 


In  the  formula,  m  is  the  mass  flow  of  the  gas,  is  the  viscosity, 
4  p  is  the  inner  and  outer  pressure  difference  of  the  sealed 
cover,  A  y  is  the  thickness  of  the  porous  metal,  S  is  its  area, 
and  Kp  is  the  coefficient  determined  by  the  porosity  of  the 
porous  metal.  If  it  is  still  required  that  the  ratio  of  the 
sealed  covers  ambient  pressure  and  the  minimum  pressure  in  the 
cover  be  0.528,  and  the  seepage  flow  is  taken  as  1.5  grams/ 
second,  then  equation  (14)  can  be  changed  to 


A*~.Pii”0.472/),J  *  (15) 

when  pq1=0.833  atmospheric  pressure,  ambient  pressure  p'=0.44 
atmospheric  pressure,  and  when  the  ambient  pressure  of  the  sealed 
cover  rises  to  an  atmospheric  pressure  the  pressure  in  the  cover 
P02=1.2601  atmospheric  pressure;  when  Pq^=0.568  atmospheric 
pressure,  pQ2=1.133  atmospheric  pressure.  Naturally,  we  can  also 
select  an  appropriate  porous  material  so  that  the  maximum 
pressure  in  the  cover  is  somewhat  lower.  For  example,  if 
p ' /Pqi=0 • 8 ,  then  when  Pg^=0.833  atmospheric  pressure,  p02=1.1236 

atmospheric  pressure;  when  pQ1=0.563  atmospheric  pressure, 
Pq2=1-°61  atmospheric  pressure. 

_4 

Under  our  conditions,  we  can  take  M=l. 85x10  grams/ 

—X  0  2 

centimeter • seconds  and  if  Kp=0. 61x10  meters  and  A  y=l 
centimeter,  then  we  take  the  porous  metal's  sectional  area 
S=1.2146  (centimeters) ^ . 
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Selection  of  the  K  value  is  based  on  the  equivalent  formula 

P 

of  the  seepage  flow  and  fine  porous  flow  [7] . 


djN, 

128 


(16) 


In  the  formula,  d.  is  the  diameter  of  the  fine  orifice  and  N.  is 

l  l 

the  number  of  fine  orifices  on  the  unit  area.  If  there  are 
twenty-five  0.1  millimeter  diameter  fine  orifices  per  square 
millimeter,  then 

K  =0 . 6136xl0~^  meters^ 

P 

Under  these  conditions,  the  aperture  rate  is  19.63  percent 
and  the  aperture  area  is  3.067  times  a  single  throttling  orifice. 

-10  2 

Therefore,  the  selection  and  placing  of  K  =0.61x10  meter 

P 

2 

porous  metal  1  centimeter  thick  and  1.2146  centimeter  sectional 
area  on  the  sealing  cover  's  shell  can  cause  pressure  in  the 
cover  to  be  maintained  in  the  li  0.3  atmospheric  pressure  range. 
Naturally,  the  inner  pressure  of  the  sealing  cover  must  be  lower 
in  an  absolute  vacuum. 

V.  Discussion  and  Conclusion 

If  the  environmental  pressure  ae-ostatic  bearing  such  as  the 
one  shown  in  fig.  1  continually  rises,  self-excited  resonance 
cannot  be  avoided  under  conditions  wherein  bearing  capacity  does 
not  influence  design.  The  resonance  source  at  this  time  is  a 
supersonic  flow  on  the  back  edge  of  the  throttling  step  on  the 
thrust  surface. 

Use  of  the  orifice  throttling  method  can  keep  the  bearing  in 
an  environment  wherein  the  sealing  cover’s  ambient  pressure 
continually  drops  and  this  will  not  cause  self-excited  resonance. 
However,  supersonic  flow  resonance  is  produced  by  the  throttling 
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orifices  on  the  sealing  cover's  shell.  According  to  experiences, 
the  influence  of  this  resonance  on  instrument  precision  can  be 
overlooked. 

Under  normal  ambient  pressure,  relatively  high  air  feed 
pressure  aerostatic  bearings  can  also  produce  self-excited 
resonance.  The  method  suggested  in  this  paper  can  also  be  used 
to  resolve  this  problem. 

If  after  test  measurements  on  the  precision  of  an  instrument 
system  it  is  necessary  to  eliminate  the  resonance  produced  by 
the  throttling  orifices  on  the  sealing  cover's  shell,  we  can 
attach  a  muffler  or  vibration  damper. 

Given  orifice  throttling,  because  the  boundary  layer  shrinks 
the  orifices  of  the  effective  sectional  area,  when  the  sealing 
cover's  ambient  pressure  is  an  atmospheric  pressure,  pressure 
pQ2  in  the  cover  must  be  slightly  higher  than  the  calculated 
results . 

To  realize  this  goal  of  damping  or  eliminating  vibration  we 
can  also  use  fine  porous  throttling  or  seeping  flow  throttling. 

As  for  fine  porous  throttling,  the  total  sectional  area  of  the 
fine  orifices  should  be  larger  than  the  sectional  area  of  a  single 
throttling  orifice.  As  for  seepage  flow  throttling,  when  the 
surrounding  ambient  pressure  is  an  atmospheric  pressure,  the 
pressure  is  an  atmospheric  pressure,  the  pressure  in  the  cover 
should  be  higher  than  orifice  throttling  under  the  same 
conditions.  If  we  select  a  porous  metal  with  K  =0 . 61xl0-^meters^ , 

2 

1  centimeter  thickness  and  a  1.2146  centimeter  area,  then  lowest 
pressure  in  the  cover  is  lower  than  the  calculated  value  when 
in  an  absolute  vacuum. 

To  sum  up,  the  use  of  orifice  throttling,  fine  porous 
throttling  and  seepage  flow  throttling  can  avoid  the  self-excited 
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resonance  of  bearings.  The  use  of  orifice  throttling  can  shift 
this  type  of  vital  self-excited  resonance  from  the  bearing  to 
the  sealing  cover's  shell  and  in  this  way  guarantee  the  normal 
operation  of  precision  instruments.  Although  there  is  no 
resonance  on  the  shell  when  seepage  flow  throttling  is  used,  yet 
the  pressure  change  range  in  the  cover  is  large.  Fine  porous 
throttling  is  a  compromise  plan  between  the  two.  However,  the 
size  of  the  orifice  apertures  and  the  number  of  fine  orifices 
must  be  determined  by  experimentation.  Because  of  throttling, 
the  bearings  always  operate  in  an  almost  atmospheric  pressure 
ambient  pressure.  This  is  not  only  advantageous  to  raising  the 
precision  of  instruments  but  can  also  greatly  raise  the 
percentage  yield  of  qualified  aerostatic  bearings. 
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A  NEW  COMPLEX  ELECTROHYDRAULIC  SERVOSYSTEM 


by  You  Jinlu 

(China  Precision  Engineering  Institute  for  Aircraft  Industry) 
Abstract 

The  servosystem  of  a  flight  motion  simulator  requires  its 

servovalve-motor  combination  to  possess  a  great  ratio  of 

continuous  adjustable  rate  (e.g.  more  than  1,100).  In  general, 

only  the  electrohydraulic  servosystem  is  available  to  the 

simulator  with  a  great  load  inertia  (e.g.  30  kilograms* 

2 

centimeters • seconds  ).  This  paper  presents  a  complex  electro¬ 
hydraulic  system,  which  increases  the  ratio  of  adjustable  rates 
from  180  to  about  1800,  and  improves  other  technical  performances 
so  as  to  meet  the  requirements  of  the  servosystem  in  the  flight 
motion  s imul a t or . 

This  paper  lists  some  ratios  of  adjustable  rates  of  domestic 
servomotors  used  in  simulators,  analyzes  the  reasons  of  the  .( ■: 
ratio  of  adjustable  rates  of  the  piston  type  motor  and  puu.s 
forward  some  measures  to  improve  it.  The  paper  also  describes 
the  principles  of  the  complex  electrohydraulic  servosystem, 
gives  the  theoretical  analysis  and  derives  the  block  diagrams 
and  transfer  functions.  Then,  it  gives  t>*»  dynamic  analysis  of 
the  complex  system  and  makes  a  comparison  with  the  original  one. 
Finally,  it  provides  the  typical  data  and  plots  of  experiments. 

I.  Preface 

Aside  from  basic  technical  requirements,  for  the  following 
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system  of  an  eiectrohydraulic  flight  motion  simulator  it  is 
especially  necessary  that  the  hydraulic  servomotor  have  a  large 
(greater  than  1,000)  continuous  adjustable  rate  and  small  low 
speed  pulsation.  Domestic  flight  motion  simulators  all  use  the 
"traditional"  eiectrohydraulic  servosystem  (called  the  original 
system) .  This  type  of  system  can  only  rely  on  a  motor  developed 
with  a  great  ratio  of  adjustable  rate  and  a  (flow)  valve  to 
match  it.  After  tests,  we  selected  a  servovalve-motor 
combination  to  accord  with  the  above  requirements.  This  design 
presented  even  greater  technical  and  manufacturing  difficulties 
for  the  motor  and  valve  which  were  not  easily  realized.  Table  1 
is  the  ratio  of  adjustable  rate  for  a  piston  type  motor  developed 
over  many  years.  Because  they  are  all  large,  the  technical 
performances  of  the  system  could  frequently  not  be  realized.  The 
eiectrohydraulic  complex  system  (called  complex  system)  used  the 
YCM-27-1  motor  and  DYSF-3Q-1  value  (they  need  not  match)  which 
were  unable  to  satisfy  the  requirements  of  the  original  system, 
and  was  thus  able  to  increase  the  servovalve-motor  combination's 
ratio  of  adjustable  rates  from  110-180  to  over  1,800.  Moreover, 
dynamic  characteristics  did  not  change  and  other  performances 
also  improved  thus  satisfying  the  requirements  of  the  system. 


(1)*  *  a  * 

(2)  *»«■.  **/«• 

C3J  "  *  * 

YMD 

20—3000 

150 

rM-00 

(4)  U3S321-3000 

100—230 

YCM-Z7-1 

(5)  <«  £101-1100 

110—1*0 

ZM-40 

10—2500 

250 

ZM-2S 

4—2500 

!  750 

Table  1 


Key:  1.  Motor  model 

2.  Rotational  speed,  rotations/minute 

3.  Ratio  of  adjustable  rate 

4.  (13  to  22) 

5.  (6  to  10) 
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II.  Means  For  Improving  the  Piston  Type  Motor 

1 .  The  Friction  Torques  of  the  Running  Components  in  a 
Piston  Type  Motor  is  the  Major  Factor  Influencing 
the  Motor ' s  Low  Speed  Performance 

The  total  friction  torque  77  [1]  of  the  running  components 
in  the  motor  is 

T,ocDm{p^pt)  a ) 

and  the  running  friction  torque  T^  [1]  is 

T, - ij-QZUP.  +  P.)  (2) 

t 

In  the  formula,  0  ,  D  ,  P,  and  P_  are  the  motor's  rotating  speed 

m  m  1  2 

displacement  and  pressures  of  the  two  chambers?  in 

P.+P_=P  ,P„  is  the  power  source  pressure?  and  C.  and  C_  are  the 
J-  2  s  s  fs  f 

static  and  running  friction  coefficients  in  the  motor.  Cf  and  Cfg 

t 

change  in  accordance  with  0  (see  fig.  1) .  The  difference 
between  the  static  running  friction  torques  is  the  reason  for  the 
production  of  the  motor's  (very  low)  starting  rotating  speed. 
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Fig.  1  Starting  and  Running  Friction  Torques  Due  to  Pressure 
in  the  Motor  Chambers 

Key:  1.  Motor's  rotating  speed  (rotations/ 
second 

We  can  see  from  fig.  1  and  formulas  (1)  and  (2)  that  as  P2 

increases  the  starting  rotating  speed  of  the  motor  rises. 

Table  2  is  typical  test  data  for  the  starting  rotating  speed  of 

the  YCM-27-1  motor.  We  can  obtain  the  same  conclusions  from 

table  2  and  reference  [1],  4-4.  At  low  speed  P-^P-j  is  verY  small 

and  thus  lowering  P  is  one  effective  method  for  improving  the 

s 

motor's  low  speed. 


**/**• 

|  *#/■**! 

4)  **«*.  */*(5) 

so 

T.S 

l.S 

3/65 

(0 

10 

1.5 

6/65 

100 

It 

3 

7/65 

no 

H 

51 

« 

0/57 

Table  2 


Table  2 


2 

Key:  1.  P  ,  kilograms/centimeter 
s  2 

2.  P  ,  kilograms/centimeter 

1  2 

3.  P  ,  kilograms/centimeter 

2  2 

4.  P^-P^/  kilograms/centimeter 

5.  Starting  rotating  speed,  rotations/ 

second 


2.  The  Motor's  Lowered  Energy  Source  Pressure  During  Low 
Speed  Does  Not  Influence  the  Motor 1 s  Output  Torque 
(Specified  Inertia  Time  of  Load  Inertia) 


The  torque  T  of  the  external  load  conveyed  by  the  motor  is 
L 

[i] 


TL~DmPL-CtDj^^~  C,Z).(P,  +Pt)  ~jqT.  (  3  ) 


In  the  formula,  PL=P1-P2;  Ca  is  the  damping  coefficient,  and  tc 

is  the  absolute  viscosity  of  the  liquid.  The  mean  flow  <?L  of 
the  motor  is  [  1] 


yL  (  4  ) 


C.  and  C  are  the  internal  and  external  leakacre  coefficients 
lm  em 

of  the  motor.  At  the  same  time,  is 


Qu-C<trx^±  (p,-po  (5) 


W  and  are  the  valve  opening  area  gradient  and  opening  travel? 

and  P  is  the  density  of  the  liquid.  By  solving  formulas  (4)  and 

(5)  we  obtain  PT  letting  K=C  .WX /p”  ;K  =2C  D  ?  C  =C.  +C  f2, 

L  a  v  i  m  m  m  im  em' 
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afterwards  we  substituted  P  value  into  formula  (3)  and  obtained 

Li 

formula  (6) 


CfD.CPt+PJ-^s-T. 


(•) 


At  low  speed,  when  the  valve-motor  combination  is  under  the 

same  inertial  load,  &  is  identical  and  P  decreases.  The  first 

m  s  F  N  1 

item  on  the  right  side  of  formula  (6)  decreases  very  little  '  , 

the  third  item  forms  a  ratio  with  P,+P  =P  which  decreases  and 

12S 

during  low  speed  P  ^  PT  .  Therefore,  Tt  [fornula  (6)]  does  not 

S  Li  Li 

decrease  because  of  the  decrease  of  P  and  T_  (with  the  same 

S  L 

inertia)  is  only  related  to  &  and  is  not  influenced  by  P 

m  s 

(naturally,  there  should  be  sufficient  P  at  the  time  to 

s  , 

guarantee  the  motor's  displacement  necessary  at  this  0  ). 

m 

Table  3  lists  experimental  data.  Therefore,  if  we  use  the  P  which 

rises  and  falls  with  the  high  and  low  of  0  ,  that  is  to  sav,  use 

m 

the  Pg  which  rises  and  falls  with  the  size  of  the  input  speed 
signal  uin«  this  can  improve  the  low  speed  performance  of  the 
valve-motor  combination. 


F.N.l  Actually  in  P^  4  —*Q.  4- — *  ,  in  order  for  0  not  to 

s  L  n\  m 

change,  Xv  must  be  enlarged.  However,  *  Pg  drops  a  little 

and  C  «  1.  Therefore,  the  first  item  on  the  right  side  of 

formuTa  (6)  decreases  even  less  with  the  P  drop. 

s 
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(Up,**/**1 

20 

20 

42 

60 

70 

90 

100 

110 

120 

130 

(2)e*»/» 

1/4 

1/7 

1/9 

1/7 

1/7 

1/8 

1/8 

1/7.5 

1/S 

1/6 

(3)F«=FS 

8.2 

7.0 

9 

10 

8.6 

10 

9 

t 

9.2 

9.2 

(4  */**» 

IS 

20 

30 

41 

SO 

60 

70 

1/24 

1/20 

1/27 

1/21 

1/22 

1/19 

1/17 

(6)FOT* 

3.4 

3.3 

3.4 

2.6 

2 

.  5 

2 

(71«  *««**««.  RUfi*TLnj]F **. 


Table  3  Experimental  Data  of  Function  T^ (P  )  (Load  and  @m 
Unchanged) 

2 

Key:  1.  Pg  kilograms/centimeter 

2 .  &  rotations/second 

m 

3.  F  (F.N.l)  kilograms  2 

4 .  P  kilograms/centimeter 

s 

5.  9  rotations/second 

m 

6.  F  (F.N.l)  kilograms 

7.  F.N.l  The  arm  was  the  same  during  testing 
and  therefore  torque  Tt  was  indicated  by 
force  F 


3.  Analysis  of  the  Valve-Motor  Combination  With  P  Rising 
and  Falling  In  Accordance  with  the  High  and  Low  of  ^m 


The  flow  (?.  of  the  ideal  zero  (square)  four  opening  flow 
Li 

valve  is  [  1] 


For  convenience  of  deduction,  if  ■>  0,  therefore 

a.-rjrx./t  •  ( 5 ) 

As  shovm  in  fig.  2(a),  we  let 

P,=Po  =  lJ„+\P.  <7> 

In  the  formula,  P  and  A  P  are  the  constant  and  variable 

so  s 

parts  of  the  variable  pressure. 


(2)  <*» 
w 


Fig..  2  Schematic  Diagram  of  Flow  Servovalve-Motor  With  Pressure 
Control  Servo 

Key:  1.  Control  current  I 
2.  Pressure  valve  p 


When  formula  (7)  is  substituted  into  formula  (5) ,  we  obtain 

QU-CJVX.J- 1  (P,„  +  \P,  -  PL)  ~K'„  (  8  ) 

In  the  formula,  K'  =C,wV  P  and  =pT  -  4  We  can  see 

qo  a  s  Li  Li  s 

that  the  flow  valve  can  use  the  1^  equivalent  for  its  flow 

because  the  P^  is  limited.  Moreover,  K^o^C^w  a/  Pg  p  shrinks  to 
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k'  and  P  changes  to  P  as  shown  in  fig.  2(b).  Because 
qo  s  so 

PJ  ^  P  ,  under  ideal  conditons  P'  is  close  to  a  constant.  When 

Jb  ^  Xj  Li 

we  compare  formulas  (5)  and  (8)  we  know  that  function 
Q.  (X  )  [or  Q.  (I*)]  is  even  closer  to  being  linear  and 

Li  V  Li 

therefore  the  pressure  control  flow  valve  is  even  more  suitable 
as  a  control  element. 

We  can  see  from  formula  (8)  that  the  analysis  of  the  valve- 

motor  combination  in  reference  [1]  is  suitable  for  use  here. 

However,  it  must  use  P  .  PJ  ,  Pi,  PJ  and  K '  to  replace  P  P. 

so  l  2  L  qo  s ,  i , 

P_,  Pr  and  K  .  We  can  also  obtain  from  its  dynamic  analysis: 

2  b  ao 

(1)  because  -VI -P  J /  P  is  close  to  a  constant,  this  causes  the 

'  L  SO 

fixed  position  spring  effect  of  the  steady  flow  force  in  the 
flow  valve  and  the  non-linearity  of  the  damping  coefficient 
produced  by  the  transient  flow  force  to  decrease;  (2)  the  liquid 
flow  force  produced  by  the  liquid  flow  impact  baffler  on  the 
loaded  spring  effect  produced  on  the  baffler  decreased  which  was 
advantageous  to  the  stability  of  the  baffle  valve;  (3)  the  static 
state  rigidity  of  the  valve-motor  combination  decreased  little 
with  the  increase  of  ?L  [  1]  . 

III.  The  Complex  Electrohydraulic  Servosystem 

1.  Function  Principles  of  the  Complex  Electrohydraulic 
System 

See  fig.  3  for  the  functioaal  diagram  of  the  complex 
electrohydraulic  (velocity)  servosystem. 


r 


Fig.  3  Functional  Diagram  of  the  Complex  Electrohydraulic 
(Velocity)  Servosystem 

Key:  1.  Amplification 

2.  Operation  selection 

3.  Correction 

4 .  Synchronous 

5.  Correction 

6.  Amplification 


The  complex  system  allows  the  valve-motor  combination  to 
operate  under  variable  oil  pressure  Pg  .  Pfi  is  realized  by  the 
pressure  control  servovalve  (called  the  pressure  valve) .  P^ 
forms  a  direct  ratio  with  control  current  Ip  [2,3].  The  frequency 
band  width  of  the  pressure  valve  is  about  100  cycles/second  and 
the  influence  of  its  dynamic  characteristics  on  the  flight 
motion  system  with  a  frequency  band  width  of  10  cycles/second 
can  be  overlooked.  The  complex  system  synchronizes  the  control 
pressure  valve  and  flow  valve  by  the  given  signal  u^n  (see  fig. 3). 
In  this  way,  when  there  is  a  small  signal  the  flow  valve  operates 
in  low  pressure,  and  with  a  small  opening  (small  flow)  the  motor 
operates  at  a  low  speed;  when  there  is  a  large  signal  the  flow 
valve  operates  in  a  high  oil  pressure  and  with  a  large  opening 
(large  flow)  the  motor  operates  at  a  high  speeds  Thus,  the  low 
speed  performance  and  dynamic  performance  of  the  valve-motor 
combination  satisfies  the  requirements  of  the  system.  The 
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stability  performance  of  the  system  is  the  same  as  that  of  the 
original  system  because  the  open-loop  segment's  stability  does 
not  have  any  influence  on  the  closed-looped  loop  included  in  the 
complex  system  [4] . 

2.  Block  Diagram  and  Transfer  Functions  of  Complex 
Servosystem  (Velocity  System) 

(1)  Block  Diagram  of  Flow  Valve- (Valvecontrolled) 

Motor  With  Pressure  Control  Servo 


We  can  infer  the  block  diagram  of  flow  valve- (valvecontrolled) 
motor  with  pressure  control  servo  from  formula  (8)  and  motor 
torque  equation  [1].  PT-^P  replaces  P'  and  after  structure 

transformation  we  obtain  fig.  4. 


Fig.  4  Block  Diagram  of  (Flow-)  Valve  (Valve-Controlled)  Motor 
With  Pressure  Control  Servo 

In  the  figure,  G q  is  the  flow  valve's  transfer  function,  the 
valve  frequency  band  is  within  100-200  cycles/second,  for  this 
system  C&  1/(T q  S  +1)  ,  and  is  a  time  constant.  The  other 

symbols  in  the  figure  are  the  basic  parameters  of  the  valve  and 

motor  [1].  Letting  K  =K  +C.  +V.S/4/6’  ,  we  can  see  that  because 

a  c  im  t  e 

the  pressure  control's  controlling  effect  on  the  flow  valve's 
flow  is  equivalent  the  control  current  1^  is  ^  Ps^KqoG<?  * 

From  the  above  mentioned  ^ P  =G  I  , G  =K  / (T  S+l) .  G  is  the 

s  p  p.  p  p  p  p 

pressure  valve's  transfer  function;  is  the  pressure  gain; 

and  T  is  the  time  constant.  When  ^P  =1  G  is  substituted  into 
p  s  p  p 


the  I ,  formula  we  obtain  the  equivalent  transfer  function  G, 
a  a 

from  I  to  I, 
p  d 

~  -  (K*+Ctm  +  V,S/i&)K.(T0S+i) 

Gt  I,  KtjTJ+l) 

Thus,  we  obtain  the  block  diagram  of  the  complex  (velocity) 
system  (see  fig.  5) . 


U,n 


Fig.  5  Block  Diagram  of  Complex  (Velocity)  Servosystem 

In  the  figure,  is  the  transfer  function  of  the  corrected 

segment,  and  H^=H^ (T^S+1 ) .  is  the  amplification  coefficient? 

2  2 

T.  is  the  time  constant;  G  =K  /(t„  S  +  2  iT„S+l)  is  the  motor's 
1  M  M  ’  M  ^  M 

transfer  function  and  K^-l/D  ;  KR  is  the  amplification 

coefficient  of  CYD-2.7. 

(2)  Transfer  Function  of  Electrohydraulic  (Velocity) 

Closed  System 

After  the  complex  system's  closed-loop  transfer  function  [5] 
written  from  fig.  5  is  revised,  it  is 


(9) 


!  KuJZK.  +  Cim  +  V.S/&rXT9S+l)(riS+l)H{K,+  K.KUTtS+i)'i 

CK.Ki.KuKHJ  +  .\(T9S  + 1  )(r&&  +  2iTuS-r  l) )  (1\S  + ij 


Key:  1.  Complex 

The  closed-loop  transfer  function  of  the  original 

system  (velocity  system)  is  (see  the  inner  part  of  the  broken 
line  of  the  block  diagram  in  fig.  5) 


.  K.K^G0GJ/N 
**^\+K.K~G9GuJ/1 


K.  K.  K  I 

KiKhKuKhJ  +N('T9S-i- 1) (TiiS1  -+- 2£7'u.S+ 1) 


Key:  1.  Original 


3.  Analysis  of  the  Dynamic  Characteristics  of  an 

Electrchydraulic  Complex  (Velocity)  Closed  System 


The  characteristic  equation  of  the  complex  closed  system 
obtained  from  formula  (9)  is 

iKlK’„KuKHJ  +  N(TQS+l)Crj,St  +  2lTllS+l'):(T'S+\)- 0  (11) 

The  characteristic  equation  of  the  original  system's  (velocity 
system)  closed  system  is  obtained  from  equation  (10) 

KtK^iKMJ  +  NtToS+l-HTlS'  +  ZlTMS+l)*  0  (12) 

Equations  (11)  and  (12)  separately  show  the  dynamic 
characteristics  of  two  types  of  closed  systems  and  there  is 
only  an  extreme  point  s=  -  ~  difference  between  them.  From 
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on 


the  point  of  view  of  the  root  locus  [5,6] ,  because  ^  on 

the  S  complex  plane,  when  compared  to  the  pair  of  conjugate 

extreme  points  of  the  motor,  extreme  point  s=  -  is  far  from 

TP 

the  imaginary  axis  and  its  influence  on  the  system's  dynamic 

characteristics  can  be  overlooked.  The  dynamic  characteristics 

of  these  two  types  of  systems  are  mainly  derived  from  the 

conjugate  extreme  points  of  the  motor,  in  other  words,  the 

dynamic  characteristics  of  the  two  types  of  systems  are 

identical.  We  know  from  eauation  (9)  that  after  H ' ,  K  and  K' 

1  p  qo 

are  suitably  disposed,  the  <^cnmpiex  numerator  arid  the  third 

degree  factor  of  the  denominator  cancel  eachother.  In  this  way, 
the  dynamic  characteristics  of  the  complex  system  are  improved 
even  more . 


IV.  Comparison  Test  of  the  Complex  Electrohydraulic  (Velocity) 
System  and  the  Original  (Velocity)  System  (Equivalent  Load 
of  System  Band  and  Flight  Motion  Simulator  Load  Equal 
Inertia) 


Two  rounds  of  multiple  open-loop  and  closed-loop  static 
tests  and  dynamic  tests  were  carried  out  on  the  two  systems 
composed  of  different  DYSF-3 Q  -IP  and  DYSF-3P  (flow)  valves  and 
the  YCM-27 -1  motor  (after  systems  tests  this  motor  was  found  to 
be  unable  to  satisfy  the  technical  index  of  the  original  system) 
with  DYSF-1P  and  DYSF-3P  pressure  valves  and  CYD-2.7  and  ZCF-16 
velocity  (they  need  not  match) .  Table  4  is  typical  data  of  two 
rounds  of  open-loop  low  speed  tests  on  the  two  types  of  systems. 


Table  4 


Key:  1.  Set 

2.  Original  system 

3.  Complex  system 

4.  Lowest  rotatinq  speed 

5 .  Rotations/minute 


We  can  see  that  the  complex  system's  valve-motor  combination 
ratio  of  adjustable  rate  enlarged  more  than  ten  fold  and 
satisfied  the  requirements  of  the  flight  motion  simulator.  Figs 
6  and  7  are  two  rounds  of  open-loop  dynamic  tests  of  the 
amplitude  (phase) -frequency  responses  of  small  signals  in  the 
two  types  of  systems  (the  signals  are  rated  values  of  0.4—16 
percent  and  1.6—68  percent).  In  the  second  round  of  tests,  when 
the  frequency  is  0.1  cycles/second,  the  open-loop  amplification 
multiple  is:  the  complex  system  is  130  1/second  ,  and  the 
original  system  only  reaches  89  1/second  because  of  instability 


Fig.  6  Amplitude  (Phase) -Frequency  Responses  of  Small  Signals 
in  the  Open-Loop  System 
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Flow  valve  DYSF-3 $ -1 (08*^;  pressure  valve  DYSF-1P; 

/  M\ 

YCM-27-1  (03  ) motor?  16  percent  of  amplitude's  rated 

value,  K  „  1/50  1/second;  K  __  1/50  1/second 
on  otr 

Key:  1.  Amplitude  (decibels) 

2.  Frequency  (cycles/second) 

3.  Phase  (degrees) 

4.  Amplitude-frequency  responses  of 

original  system 

5.  Amplitude  frequency  responses  of 

complex  system 

6.  Phase-frequency  responses  of  original 

system 

7.  Phase-frequency  responses  of  complex 

system 


(2  )**<*/*> 


Fig.  7  Amplitude  (Phase) -Frequency  Responses  of  Small  Signals 
in  the  Open-Loop  Systems 

Flow  valve  DYSF-3  <^-l  ^04**);  pressure  valve  DYSF-3P 

YCM-27-1  (05*)  motor;  7  percent  of  amplitude’s  rated 
value 
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Key:  1 
2 

3 

4 

5 

6 
7 

Fig.  8  is  two  rounds  o 
two  types  of  closed-lo 
between  each  wave  form 
unrelated) . 

(1)  l.B** 

( 2 )  s*  « 


Fig.  8  Output  Signal  Wave  Forms  of  (Velocity)  Closed -Loop 
Systems 

Key:  1.  Original  system  u^n=  3  volts,  Uin  is 

the  sine  value 

2.  5  cycles/second 

3.  15  cycles/second 

4.  20  cycles/second 

5.  Complex  system,  U.  =2  volts,  U.  is 

in  in 

the  sine  value 

6.  5  cycles/second 

7.  30  cycles/second 

8.  35  cycles/second 


Fig.  9  is  the  amplitude  (phase) -frequency  response  of  big  signals 
in  the  open-loop  complex  system  (68  percent  of  the  rated  value) . 
The  frequency  band  width  drop  from  the  35  cycles/second  of  the 
small  signals  to  30  cycles/second  is  rational  [7] . 


^2)  »(W») 


Fig.  9  Amplitude  (Phase_-Frequencv  Response  of  Big  Signals 
in  the  Open-Loop  Complex  System 

Flow  valve  DYSF-3(^ -1  (04^)  ;  pressure  valve  DYSF-3P; 

YCM-27-1  f05^)  motor;  16  percent  of  amplitude's 
rated  value. 
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Key:  1.  Amplitude  (decibels) 

2.  Frequency  (cycles/second) 

3.  Phase  (degrees) 

4.  Amplitude-frequency  responses 

5.  Phase-frequency  responses 

Figs.  6  and  7  show  that  the  dynamic  characteristics  of  the  two 
systems  are  the  same  and  that  the  damping  of  the  complex 
system  is  larger.  Fig.  8  shows  that  the  complex  system's  output 
distortion  is  small,  its  dead  area  is  small  and  its  damping 
coefficient  is  large  so  that  controllability  and  stability  are 
both  good. 

To  sum  up,  when  similar  motors  and  valves  are  used,  the 
complex  system,  as  compared  to  the  original  system,  can  greatly 
improve  the  low  speed  performance  of  the  valve-motor 
combination  and  cause  the  ratio  of  the  adjustable  rate  of  the 
valve-motor  to  enlarge  more  than  ten  fold,  the  system  to  be 
non-linear,  the  dead  area  to  be  relatively  small  and  the  system's 
damping  coefficient  to  be  relatively  large,  when  the  open-loop 
amplification  multiple  is  the  same,  the  system's  stability  is 
good  and  thus  it  is  easy  to  realize  the  technical  performance  of 
the  flight  motion  simulator's  servosystem.  This  decreases  the 
harsh  technical  requirements  for  the  motor  and  valve  and 
decreases  the  difficulty  of  developing  motors.  It  also  raises  the 
technical  performance  of  the  system. 
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A  COMPUTATIONAL  GEOMETRY  METHOD  FOR  BLADE  SPACE  GEOMETRIC 
DESIGN 
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Zhao  Yuqi,  Zhang  Tingxiong  and  Xiao  Hongen 
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Abstract 

This  paper  presents  a  new  method  of  computer-aided  blade 
geometric  design  which  combines  several  typical  methods  of 
computational  geometry  developed  recently  (including  Bezier 
patches,  B-spline  patches  and  Bezier-B-spline  mixed  patches 
etc.)  with  the  space  shaping  for  aerojet  blades.  An  integral 
mathematical  model  for  a  blade  configuration  is  obtained  and 
the  curved  blade  is  jointed  by  a  set  of  3  times  2  order  mixed 
surface  patches.  Thus,  it  is  advantageous  to  the  computational 
analysis  of  the  blade's  three  dimensional  flow  field  performance 
and  limited  strength  as  well  as  combining  with  numerical  control 
processing.  Actual  cases  show  that  this  method  is  simple,  has  a 
small  amount  of  computations  and  is  effective  for  the  design  of 
blades  with  large  curvatures. 

I .  Preface 

Blade  shape  design  is  an  important  key  in  computer-aided 
blade  design  and  is  closely  connected  to  the  blade’s 
aerodynamic  and  strength  computations.  In  the  past,  plane  surface 
shaping  was  commonly  used  in  engineering  and  later  they  carried 
out  the  method  of  space  superposition  for  the  given  profilers 
center  of  gravity  law  along  the  high  point  of  the  blade. NThere 
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is  no  integral  mathematical  model  for  the  blade  and  processing 
is  complex  which  presents  difficulties  for  computer-aided 
design  and  manufacture. 


Present  computational  geometry  has  attained  remarkable 
advancements  in  configuration  designs [  1-4] .  This  paper  attempts 
to  find  and  investigate  the  application  of  certain  of  the  newest 
achievements  [4,5]  in  the  configuration  design  of  aerojet  blades 
as  well  as  the  manufacturing  laws  for  high  load  large  curved 
turbine  blades.  A  certain  number  of  design  examples  have  shown 
that  the  space  shaping  method  for  blades  suggested  in  this  paper 
is  simple,  flexible,  directly  perceived  and  requires  a  small 
amount  of  computations.  The  blade  profile  has  an  integral 
mathematical  formula  which  can  conveniently  find  the  contour  of 
an  arbitrary  section,  is  convenient  for  computing  the  aerodynamics 
of  the  blade's  three-dimensional  flow  and  analyzing  the  finite 
element's  strength.  Moreover,  processing  is  simple  and  it  is  not 
necessary  to  carry  out  blade  fitting  [6] .  This  method  is  also 
easily  extended  to  conical  and  flow  surfaces. 

II.  Mathematical  Model 

Although  the  curvature  of  the  blade  is  large  chordwise,  yet 
the  blade's  high  direction  becomes  level  approaching  a  straight 
line.  Therefore,  the  twisted  space  patch  is  close  to  a  ruled 
surface  or  is  a  ruled  surface.  For  this  reason,  we  use  a  surface 
patch  along  the  high  part  of  the  blade  and  use  any  surface  patch 
chordwise  to  describe  the  blade  profile.  Moreover,  in  the 
surface  patch  used,  the  u  line  uses  the  parametric  cubic  curve 
and  the  v  line  uses  the  parametric  quadratic  curve,  that  is,  the 
3  times  2  order  mixed  surface  patch.  They  can  utilize  a  unified 
mathematical  formula: 
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In  the  formula,  [B]  is  the  characteristic  polyhedron  vertex 
position  vector  matrix;  i  is  surface  patch  sequence  number; 
and  k  is  the  mixed  surface  patch  type  index. 

When  k=l ,  the  u  line  of  the  surface  patch  is  a  uniform 
cubic  B-spline  curve  and  should  take 


(1»> 


When  k=2,  the  u  line  is  the  X-spline  suggested  by  Sun  Jiachang 
and  at  this  time 


-  1 


CM], 


1 

0 

0 

0 


(lb) 


When  k=3,  the  u  line  is  the  cubic  Bezier  curve  and  at  this  time 
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When  the  v  line  is  a  quadratic  Bezier  curve  and  there  are  no 
interpolation  requirements,  at  this  time 


)-[-2  I  ol 

Cl  0  0  J 


In  order  for  the  contour  formed  from  the  surface  to  satisfy 
points  E  and  H,  we  used  EF  and  HG  as  its  tangent.  Therefore,  its 
vertex  must  be  extrapolated  according  to  the  parabolic  qualities 
(see  fig.  1) . 


/iLwi 


£'mT  > 

x  J  ✓- 


Fig.  1  3  Times  2  Order  B4zier  Patches 


After  deduction 
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Naturally,  after  this  type  of  processing,  the  v  line  of 
the  surface  is  the  quadratic  Bezier  curve  of  the  interpolation 
requirement.  Moreover,  the  u  line  in  the  middle  section  must 
be  tangent  to  the  given  leading  and  trailing  edge  points  of 
tangency . 

III.  Design  Method 

How  do  we  provide  the  characteristic  polyhedron  vertex 
position  vector  and  how  many  patches  are  used  to  describe  the 
blade  prfile?  This  is  the  crux  of  this  method. 

We  divided  the  apex  position  vortex  of  the  characteristic 
polyhedron  in  space  into  three  basic  sections.  Firstly,  based 
on  past  blade  shaping  experience  [7]  we  gave  the  basic  point 
of  the  characteristic  polygon  on  the  basic  section  and 
afterwards  gave  the  necessary  control  points  based  on  the  type 
k  of  the  computational  geometry  method  and  went  through  revision 
for  final  complete  shaping.  For  this  reason,  we  used  the  two 
paths  ard  three  types  method  for  specific  shaping. 

1.  The  Bezier  Method  of  Blade  Shaping 

This  is  a  method  which  uses  a  small  amount  of  surface  patch 
to  individually  describe  the  curve  and  back  profile  of  the  blade. 
Because  the  curvature  changes  of  the  blade  curve 1 s  contour  are 
relatively  uniform,  we  used  3  times  2  Bezier  surface  patches? 
the  curvature  changes  of  the  back  profile  of  the  blade  were 
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relatively  large  and  therefore  was  joined  by  two  Bezier  surface 
patches  (see  fig.  2) . 


Fig.  2  A  Surface  Built  by  Two  3  Times  2  Order  Brezier  Patches 

Key :  1 .  Root 

2.  Middle 

3 .  Apex 

4.  Surface  patch  I 

5.  Surface  patch  II 


To  satisfy  joining  C  continuously,  we  used  strong  ample 
conditions,  and  not  including 


(  /  =*  0,  1 ,  2)  (2) 


the  collinear  three  points  must  also  satisfy: 


In  the  formula,  A  is  the  proportionality  constant. 


When  using  the  Bezier  method  for  blade  shaping  it  is 
necessary  to  give  the  boundary  conditions  of  the  surface  patch. 
The  two  joining  point  position  vectors  and  common  tangent 
directions  must  also  be  determined  for  the  back  of  the  blade. 
Generally,  this  joining  point  has  the  straight  line  section  in 
back  of  the  blade  as  the  starting  point  and  the  tangent  direction 
is  its  linear  direction.  If  there  is  no  linear  section  then  this 
point  is  selected  on  the  throat  tangent  point  G  area.  Its  tangent 
direction  is  determined  by  the  tail  end  corner  on  the  back  of  the 
blade.  In  short,  the  determination  of  the  cubic  Bezier  curve  on 
the  basic  section  is  similar  to  the  two  point  two  tangent  method 
commonly  seen  in  engineering  curve  shaping.  The  different  ones 
still  must  determine  the  other  two  apex  B^  and  B2j  along  the 

tangent  direction  of  the  two  extreme  points.  This  paper  used  a 
large  curvature  blade  back  and  blade  curve  as  references.  They 
were  determined  by  test  gathering  method  and  the  results  were 
excellent . 

2.  The  B-Spline  Method  of  Blade  Shaping 

The  whole  blade  profile  is  composed  of  a  set  of  3  times  2 

order  B-spline-Bezier  mixed  surface  p.  tches  and  the  joining 

2 

area  maintains  C  continuously  (see  fig.  3) . 
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Fig.  3  Shaping  by  Bezier-B-Spline  Mixed  Patches 

Key:  1.  Base 

2.  Middle 

3 .  Apex 

The  central  problem  of  this  method  is  giving  the  vertex 
position  vortices  of  the  sealed  polygons  on  the  three  basic 
sections. 

Naturally,  it  is  not  enough  to  only  rely  on  the  polygon 
formed  from  the  basic  points  to  describe  the  blade  profile.  It 
is  also  necessary  to  give  the  control  points  based  on  the 
special  features  of  the  different  blade  profiles  as  well  as  the 
specific  requirements  for  each  shaping.  Because  there  are  so 
many  blade  profile  types  there  are  differences  in  each  of  the 


shaping  requirements.  Therefore,  it  is  very  difficult  to  give 
a  uniform  pattern.  We  will  only  list  a  few  principles  below  as 
a  reference  for  concrete  implementation. 

(1)  It  is  best  that  the  characteristic  polygon  in  the  leading 
and  trailing  edge  areas  cause  the  curve  approximation  to  be 
round ; 


(2)  A  control  point  can  be  selected  in  the  large  curvature 
area  near  the  leading  edge  of  the  blade  profile; 


(3)  If  a^  is  the  length  of  the  side  of  the  characteristic 
polygon  and  a.^  is  the  corner  of  the  adjacent  side,  then 
aitg(ai/2)  can  roughly  realize  the  curve  radius  of  this  area's 
B-spline  curve.  Thus,  we  can  use  the  length  of  the  side  and  the 
corner  of  the  control  characteristic  polygon  to  satisfy  the 
changing  tendency  of  the  blade  profile's  curvature.  For  example, 
if  the  blade  profile's  curvature  radius  requires  monotone 
increasing,  then  it  is  necessary  for  the  series  formed  from 

aitg^2  )  to  also  monotone  increase. 


3.  The  X-Spline  Method  of  Blade  Shaping 

For  convenience  of  narration,  we  called  the  spline  method 
proposed  in  reference  [5]  the  "X-spline."  In  this  way,  the 
whole  blade  profile  is  composed  from  a  set  of  3  times  2  order 
X-spline-Bezier  mixed  surfaces  and  the  joining  area  maintains 

C1  continuously.  Its  special  feature  is  that  there  is 
interpolation  in  the  u  line.  It  is  only  necessary  to  properly 
select  the  apex  of  the  polygon  for  the  contour  to  be  able  to 
use  the  control  points.  The  other  situations  are  the  same  as 
the  B-spline  shaping  method  and  need  not  be  explained  here. 

Because  the  blade's  configuration  design  not  only  seeks  the 
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"beauty"  of  shaping  but  more  importantly  must  satisfy  certain 
aerodynamic  and  mechanical  performance  requirements,  the 
method  presented  in  this  paper  is  only  a  designing  tool  for  the 
blade's  space  analytical  forming.  Whether  or  not  the  designed 
blade  is  satisfactory  it  is  still  necessary  to  calculate  the 
flow  around  the  cascade  and  evaluate  by  means  of  checking 
velocity  and  pressure  distribution.  When  compared  to  other 
shaping  methods,  this  method  is  very  flexible  in  respect  to 
modifying  the  design. 

IV.  Design  Examples 

To  guarantee  the  triangle  of  velocity  for  the  aerodynamic 
design,  we  depend  on  engineering  experience  to  first  select 
initial  geometric  parameters  for  blade  shaping  on  each  basic 
section  as  shown  in  table  1 . 
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Table  1.  Initial  Geometric  Parameters  For  Blade  Shaping 

Key:  1.  Sequence  number 

2 .  Symbol 

3 .  Name 

4 .  Leading  edge  radius 
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Key: 


5.  Trailing  edge  radius 

6.  Inlet  structural  angle 

7.  Outlet  structural  angle 

8 .  Leading  edge  wedge  angle 

9.  Trailing  edge  wedge  angle 

10.  Axial  chord 

11.  Setting  angle 

12.  Throat  radius 

13.  Grid  pitch 

14.  Exit  edge  corner  on  back  of  blade 

15.  Maxium  thickness 

16.  Original  point  on  basic  coordinate  X 

17.  Original  point  on  basic  coordinate  Y 

18.  Section  on  basic  coordinate  Z 

19.  Unit 

20.  Blade  root 

21.  Middle  of  blade 

22.  Blade  apex 

23.  Millimeters 

24.  Millimeters 

25.  Degrees 

26.  Degrees 

27.  Degrees 

28 .  Degrees 

29.  Millimeters 

30.  Degrees 

31.  Millimeters 

32.  Millimeters 

33 .  Degrees 

34.  Millimeters 

35.  Millimeters 

36.  Millimeters 

37.  Millimeters 


Based  on  the  geometric  relationship  on  the  basic  section,  it  is 
very  easy  to  find  basic  points  A,  B,  C,  D,  E,  F,  G  and  Q  of 
the  control  blade  profile  as  well  as  the  tangent  directions  of 
related  points  (see  fig.  4) . 
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Fig.  4  Initial  Geometric  Parameters  and  Basic  Points  of  a 
Blade  Profile 


The  shaping  method  for  B-spline  and  X-spline  blades  also 
requires  that  an  additional  adjustable  basic  point  be  added 
between  the  A,  F  point  and  F,  G  point  on  the  back  of  the  blade 
as  well  as  one  between  the  C,  E  point  and  E,  D  point  on  the 
curvature  of  the  blade.  The  induced  adjustable  parameters 

KJ,  KI^  . KI^  and  other  related  adjustable  parameters  are 

listed  in  table  2. 


1  *  * 

*C2I* 

«  C3J  * 

(13) 

A  <2 

"(19)  1 
*»■  * 

— (W 

•7  f 

— tm- 
!  *  n 

1 

Fl 

FMAWXMB.Ztt 

(4) 

— 

0.434 

0.338 

0.34 

2 

Fh 

F  *B.Zit 

(5) 

— 

0.2*5 

0.115 

-0.17 

*r 

F  Attft* 

(6) 

*(14 

-7.5 

-23 

-52 

4 

K0 

********* 

(7) 

— 

0.5 

0.3 

0.35 

5 

KJ 

(8} 

— 

0.7 

0.7 

0.7 

1 

KIt 

(9) 

— 

0.4 

0  34 

0.4 

7 

Kh 

*%S«A*T**ft 

(10) 

— 

045 

0.58* 

0.7 

* 

Kh 

(ID 

— 

0.5 

0.5 

0.5 

• 

Kit 

(12) 

— 

0  44 

0.5 

0.522 

Table  2  Adjustable  Parameters  of  the  Control  Polygon 


Key :  1 .  Sequence  number 

2 .  Symbol 

3 .  Name 

4.  Ratio  of  F  point  coordinates  X  and 

5.  Ratio  of  F  point  coordinates  Y  and  Bx 

6.  Oblique  angle  of  F  point 

7.  Adjustable  parameters  of  linear  section 

on  back  of  blade 

8.  Supplementary  basic  point  adjustable 

parameters 

9.  Supplementary  basic  point  adjustable 

parameters 

10.  Supplementary  basic  point  adjustable 

paramenters 

11.  Supplementary  basic  point  adjustable 

parameters 

12.  Supplementary  basic  point  adjustable 

parameters 

13.  Unit 

14.  Degrees 

15.  Blade  root 

16.  Middlo  of  blade 

17.  Blade  apex 


In  order  for  the  leading  and  trailing  edges  to  be  approximately 
round,  we  formed  a  regular  polygon  by  evenly  dividing  a  central 
angle  between  two  tangent  points  (for  example,  the  A  and  C 
points  in  fig.  4)  of  the  leading  and  trailing  edges  (as  shown  in 
fig.  5) . 


Fig.  5 
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Fig.  5  Formation  of  Polygonal  Vertices  at  the  Leading  Edge 
of  a  Blade  Profile 

This  not  only  controls  the  basic  contour  of  the  blade  profile 
from  the  basic  point  but  also  guarantees  shaping  requirements. 
Moreover,  depending  on  this,  there  is  basically  formed  a  control 
polygonal  apex.  The  three  basic  sections  are  controlled  in  this 
way,  that  is,  they  form  control  polygons  of  the  blade  profile's 
space  patches.  Limited  by  space,  this  paper  only  lists  the 
initial  geometric  parameters  and  adjustable  parameters  of  the 
X-spline  method  for  blade  shaping  in  tables  1  and  2.  See  fig.  6 
for  the  design  results. 


Fig.  6  Shaping  a  Blade  by  the  X  Spline  Method 

They  are  directly  drawn  on  a  drawing  instrument  by  computer. 
Aside  from  this,  the  B-spline  and  Bezier  methods  are  separately 
used  to  design  the  shapes  of  the  other  two  blade  profiles. 

Only  the  related  results  are  drawn  in  figs.  7  and  8. 
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Fig.  7  A  Profile  Shaped  by  the  B-Spline  Method  and  Its 
Characteristic  Polygon 


Fig.  8  Shaping  a  Blade  by  the  Bezier  Method 
V.  Conclusion 

1)  After  analysis  and  sample  calculations  we  know  that  an 
effective  new  shaping  method  can  directly  and  competently  form 
the  design  for  the  space  of  various  types  of  large  curvature 
blades.  The  whole  biade  profile  has  an  integral  mathematical 
model  which  is  very  convenient  for  the  computer-aided  design 
and  manufacture  of  blades. 


2)  Each  of  the  three  types  of  shaping  methods  have  their 
special  features.  B-spline  shaping  has  high  patch  smoothness, 
revisions  are  flexible  and  its  localized  strength  gives  people 
the  feeling  that  its  shaping  curve  is  "flexible."  The  contour 
made  by  the  Bezier  method  is  "hard,"  its  smoothness  is  good  and 
its  required  surface  patches  are  small.  Yet,  to  guarantee  a 
certain  smoothness  at  the  joining  area  of  the  surface  patches 

it  is  necessary  to  consider  special  allocation  of  the  apex. 
X-spline  shaping  lies  between  the  two  and  its  one  outstanding 
point  is  that  it  has  interpolation  like  that  of  the  Bezier 
method.  Moreover,  it  can  also  automatically  guarantee  C* 
smoothness  on  the  patch's  joining  area. 

3)  This  new  shaping  method  is  beneficial  to  three  dimensional 
flow  field  calculations,  finite  element  strength  analysis  and 
numerical  control  processing.  It  is  also  a  useful  tool  in  making 
an  integral  whole  of  the  CAD/CAM  blade.  However,  this  is  only  a 
beginning  and  we  must  await  a  large  quantity  of  design  practice 
and  the  accumulation  of  experiences.  The  improvements  are  not 
sufficient  especially  when  used  for  flow  surface  shaping.  The 
problem  of  how  to  give  the  characteristic  polyhedron  apex 
position  vortex  based  on  requirements  also  requires  detailed 
discussion. 
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OVERSPEED  AND  OVERTEMPERATURE  TESTS  AND  DISSECTION  EXAMINATION 
OF  THE  TURBINE  DISCS  MADE  OF  A  IRON -NICKEL  BASE  SUPERALLOY 

by  Guo  Jianting  and  Zhang  Jinyan 

(Institute  of  Metal  Research,  Academia  Sinica) 

A  35  Ni-15  Cr  type  iron-nickel  base  superalloy  is  a  sort  of 
engine  turbine  disc  material  developed  in  China.,  It  possesses 
outstanding  comprehensive  properties  and  microstructure 
stability.  Therefore,  it  can  meet  the  requirements  of  the 
turbine  discs  of  aeroengines  for  long  time  use  below  650— 700°C. 
For  e.  new  type  of  material  to  be  formally  used  in  an  engine  it 
is  necessary  to  go  through  long  time  test  runs  and  test  flights 
However,  these  types  of  tests  are  expensive  and  long,  and  so 
to  examine  the  turbine  disc's  body  and  tenon  it  is  reasonable 
to  first  conduct  overspeed  and  overtemperature  tests  on  the 
finished  disc. 

These  tests  were  conducted  on  Chinese  made  vertical  testers. 
Six  thermoelectric  couples  are  installed  on  different  positions 
of  the  turbine  disc  to  measure  the  temperature  field.  A 
photoelectric  velocity  measurement  apparatus  and  gear  velocity 
measurement  apparatus  are  used  simultaneously  to  measure 
rotation.  To  measure  the  size  change  of  the  diameter  after  a 
turbine  disc  overspeed  test,  four  different  position  marks  were 
made  on  the  turbine  disc  prior  to  the  test.  Overtemperature 
tests  were  conducted  on  the  third  level  turbine  disc  of  a 
turboprop  engine  and  the  free  turbine  disc  of  a  turbine  axis 
engine. 

An  overspeed  test  (overspeed  15  percent)  on  the  turbine  disc 
of  the  turboprop  engine  is  first  carried  out  at  normal 
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temperature  and  after  operating  for  5  minutes  in  the  designed 
temperature  field,  the  turbine  disc  successfully  passed  the 
5  minute  15  percent  overspeed  examination.  After  examining  the 
turbine  disc  it  was  put  through  overtemperature  tests.  The 
turbine  disc  went  through  a  45°C,  5  minute  test  under  maximum 
design  temperature. 

After  the  turbine  disc  of  the  turbine  axis  engine  was  in  a 
design  temperature  field  of  10  percent  and  20  percent  overspeed, 
the  engine  was  stopped  for  cooling  examination  and  when  a 
magnifier  was  used  on  the  surface  there  were  no  abnormalities. 
Measurement  of  the  diameter  showed  an  increase  of  0.61  percent, 
overspeed  of  30  percent  and  after  being  stable  for  3  minutes, 
examinations  revealed  that  all  was  normal.  Later,  rotation 
slowly  climbed  to  over speed  of  43  percent  and  there  was  a 
"pinging"  sound.  Immediately  after  stopping  the  engine  examina¬ 
tions  showed  that  the  turbine  disc  had  broken  and  flown  out  and 
that  the  blades  on  the  disc  were  thrown  out.  Nine  mortise  of  the 
disc  were  broken  and  most  of  the  mortise  locations  had  cracks. 
There  were  relatively  light  rippled  flows  on  the  surface  of  the 
disc,  and  the  fixed  screw  holes  became  elliptical  and  most  of 
the  holes  had  cracks.  Measurement  of  the  size  of  the  disc  showed 
that  the  mean  increase  of  the  size  of  the  diamter  was  a  weak 
1.8  percent. 

We  carried  out  dissection  examination  tests  to  understand 
the  properties  of  the  overspeed  broken  disc  and  to  see  whether 
there  was  a  change  in  its  structure. 

Hardness  measurement  of  the  disc’s  longitudinal  section 
showed  that  the  hardness  value  was  basically  the  same  from  the 
center  to  the  edge.  See  fig.  1  for  the  tensile  and  shock  tests 
at  room  temperature  and  550°C. 


199 


ffjli'pBhf  - 'iff  1  "ii-  MMEMtteL 


(1)  %  5 


200  400  600 

(3)«*CD 


Fig. 


Transient  Mechanical  Properties  of  the  Broken  Disc 
Under  Over speed  Test  (Curves  -  of  Its  Forge  Disc  Broken 
Properties;  Points  -  of  the  Overspeed  Disc  Properties) 

2 

Key;  1.  (Kilograms  *meters/centimeters  ) 

2 

2.  (Kilograms/millimeters  ) 

3.  (Temperature  (°C) 


When  compared  with  a  turbine  disc  forging,  tensile  strength, 
yield  strength,  plasticity  and  toughness  are  on  the  same  level. 
See  fig.  2  for  endurance  time  at  550  and  750>°C. 


Fig. 


Stress-Rupture  Properties  of  the  Disc  Under  Overspeed 
Test  (Curves  -  of  its  Forge  Disc  Broken  Properties; 

Points  -  of  the  Over speed  Disc  Properties) 

2 

Key;  1.  Stress  (kilograms/millimeter  ) 

2.  Time  (hours) 

We  can  see  that  under  550°C,  83  kilograms/millimeter  2  conditions, 

endurance  life  does  not  drop.  Under  750°C,  30  kilogram/ 

2 

millimeter  conditions,  endurance  time  lowers. 

See  fig.  3  for  a  scanning  electron  micrograph  of  the  mortise 
fracture  of  the  turbine  disc 


Fig.  3  A  Scanning  Electron  Micrograph  of  the  Mortise  Fracture 
of  the  Over speed-Tested  Disc  (x200) 

The  fracture  has  a  crystal  sugar  type  structure  and  is  a 
typical  crystalline  fracture.  Examination  of  the  lonoitudinal 
low  multiple  structure  of  the  overspeed  broken  disc  shows  that 
the  structure  is  normal  and  there  are  no  low  multiple  alloy 
defects.  See  fig.  4  for  the  metal  phase  structure. 


Fig.  4  Metal  Phase  Structure  of  the  Overspeed -Tested  Disc 
(x500) 

Aside  from  the  Tic  and  seen  in  the  normal  structure,  there 

were  no  abnormal  phases  or  abnormal  structures.  The  quantity  of 
carbides  is  0.58  percent  which  is  the  same  as  the  quantity  of 
TiC  contained  in  the  normal  structure  of  the  35Ni-15Cr  alloy. 
There  is  about  15  percent  Y  *  phase  in  the  alloy  and  this  is  the 
major  strength  phase  of  the  alloy.  Yet,  because  its  diameter  is 
only  about  500  8,  the  metal  phase  photo  is  not  shown.  Aside  from 
this,  the  alloy  also  contains  a  minute  amount  of  C#  phase.  Long 
time  testing  of  2,000  hours  showed  that  this  quantity  of  4  phase 
had  a  very  small  effect  on  performance  and  use. 

Because  the  results  of  the  overspeed  and  overtemperature 
tests  were  excellent,  after  all  of  the  turbine  discs  were 
installed  in  the  enqine  they  underwent  multiple  trial  runs. 
Several  types  of  engines  withstood  all  of  the  long-time  trial 
flight  tests. 

We  can  derive  the  following  conclusions  from  the  above  test 
results : 

Turbine  discs  manufactured  from  35Ni-15Cr  type  iron-nickel 
base  superalloy  successfully  withstood  10,  15,  20  and  30  percent 
overspeed  tests.  The  former  also  smoothly  underwent  45°C 
overtemperature  tests  and  the  situation  was  excellent.  From 
calculations  of  overspeed  broken  rotating  velocity,  the  latter's 
safety  coefficient  reached  to  about  1.43. 

The  transient  strength,  plasticity  and  toughness  of  the 
overspeed  broken  disc,  the  endurance  time  of  550°c  and  the 
hardness  from  the  center  to  the  edge  can  still  maintain  the 
original  levels.  Low  multiple  and  metal  phase  structures  are 
normal . 
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Zhuzhou  Engine  Research  Institute,  Dongan  Machinery  Plant  and 
Hongxiangj iang  Machine  Plant  for  their  participation  in  the 
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THE  DELEGATION  OF  CSAA  TOOK  PART  IN  THE  FIRST  INTERNATIONAL 
AEROSPACE  SYMPOSIUM  IN  PARIS 

A  delegation  of  four  people  from  the  CSAA  led  by  Shen  Yuan, 
director  of  the  CSAA,  took  part  in  the  First  International 
Aerospace  Symposium  held  in  Paris  from  June  2  to  3  of  1981. 

This  symposium  was  proposed  and  sponsored  by  the  United 
States  Department  of  Commerce,  United  States  Airforce  and  United 
States  Aerospace  Bureau,  was  agreed  upon  by  the  Aerospace 
Industry  Federation  of  France,  and  was  arranged  by  the  United 
States  Association  of  Aeronautics  and  Astronautics. The  stated 
aim  of  the  conference  was  to  "promote  the  true  international 
understanding  in  the  quality  of  aeronautics  and  astronautics 
and  provide  an  opportunity  for  beneficial  exchange  among  the 
aerospace  nations  of  the  world."  When  considering  the  "leaders 
of  the  aeronautics  circles  of  each  nation,  the  Paris 
Aeronautics  Exhibition  was  the  most  magnificent  as  well  as  the 
most  pleasing  site."  Therefore,  the  Paris  exhibition  site  was 
chosen  and  during  the  two  days  prior  to  the  exhibit  the 
symposium  was  convened  there.  The  United  States  proposed  that 
the  exhibit  hall  be  the  location  of  the  first  symposium  and 
hoped  that  this  type  of  symposium  would  be  convened  once  a  year 
in  the  future. 

The  general  membership  preparatory  committee  invited  250 
representatives  from  over  20  nations  to  participate  in  the 
symposium  including  the  United  States,  France,  England,  West 
Germany,  Japan,  Italy  etc.  Prior  to  the  meeting,  the  Soviet 
Union  had  confirmed  that  they  would  participate  but  they  did 
not  attend  the  meeting.  There  were  also  reporters  from  various 
nations . 
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The  symposium  was  divided  into  four  units  and  each  unit  was 
held  for  one  half  day.  The  speakers  were  invited  and  issued 
written  materials  before  the  meeting.  During  the  two  days  a 
total  of  25  papers  were  delivered. 

The  general  subject  of  the  first  morning  were  government 
policies  and  plans  in  aeronautics  and  astronautics.  The  speakers 
were  all  technical  leaders  of  each  nation's  aerospace  science 
and  technology  leading  organization  or  large  research  units.  The 
contents  included:  1.  a  speech  by  Luofulaisi,  Deputy  of  the 
United  States  Department  of  Aeronautics  and  Astronautics,  on 
the  policies  of  that  department;  2.  a  speech  by  Professor  Telela, 
Head  of  the  European  Bureau  of  Space  Technology,  on  the  policies 
of  that  bureau;  3.  a  speech  by  Yoshiki,  Chairman  of  the  Japanese 
Commission  on  Space  Enterprises,  on  Japan's  space  activities  and 
plans;  4.  a  speech  by  General  Smith,  Chief  of  the  General  Staff 
of  the  European  Alliance,  on  the  military  use  of  aeronautics  and 
astronautics;  5.  a  speech  by  Shen  Yuan  on  certain  aspects  of 
China's  aeronautics  and  astronautics;  6.  a  speech  by  Kelukasi, 
Chairman  of  the  World  Communications  Satellite  Association,  on  the 
use  of  satellites  for  international  communications;  7.  a  speech 
by  Haweiian-Faerke,  Chief  Engineer  of  France's  Dasuo  Company,  on 
the  problems  of  the  military  use  of  aeronautics  and  international 
policies.  Shen  Yuan’s  speech  was  well  received  by  the  members  and 
was  helpful  for  foreign  nations  to  understand  the  past  and  present 
developments  of  China's  aerospace  enterprises,  related  Chinese 
policies,  and  the  general  situation  of  China's  aeronautical 
education.  The  250  copies  of  the  speech  carried  by  our  deleaation 
were  very  quickly  distributed  and  after  the  symposium  there  was 
still  a  demand  for  them.  This  speech  was  printed  in  abstract  form 
in  the  Daily  News  of  the  Aerospace  Exhibition. 

The  main  speeches  of  the  second  unit  were  speeches  by  Caisi 
Taile,  Vice  President  of  the  Being  Company,  on  the  development  of 
transporters  in  the  1980 's;  Simoerwote,  Marshal  of  the  Royal 
Airforce  and  Advisor  to  the  British  Aerospace  Company,  on  the 
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future  military  use  of  aircraft?  and  Professor  Bosangkentang, 
Director  of  the  International  Cooperation  Department  of  France's 
Aerospace  Research  Institute,  on  the  institute's  research  work. 
The  main  speeches  of  the  third  unit  were  by  Bulacai,  Director 
of  the  Applied  Plans  Department  of  France's  Aerospace  Research 
Institute,  on  the  prospects  of  space  technology  for  ground 
observations?  Professor  Dengnisi,  Chairman  of  the  Space  Research 
Committee,  on  space  science  and  international  plans?  Cailierlike, 
Vice  President  of  the  Luokeweier  International  Company,  on  space 
transport  systems.  It  was  originally  planned  that  Liuming, 
astronaut  of  the  "Salute  6",  was  to  give  a  speech  on  the  use  of 
orbiting  space  stations  and  their  effects  on  the  people's 
economy  but  this  did  not  come  about.  The  first  half  of  the  fourth 
unit  were  speeches  by  Caiqiewalei,  Vice  President  of  the  French 
Aerospace  Company,  on  civil  aviation?  Cailuojiesi,  Vice  President 
in  charge  of  international  plans  for  the  F-16  at  the  General 
Motors  Company  of  the  United  States,  on  the  prospects  of  the 
military  use  of  aircraft  etc..  The  second  half  were  speeches  by 
the  American  astronauts  Young  and  Kilibenzuo,  on  the  first 
orbital  flight  of  a  spacecraft. 

After  the  symposium,  the  Chinese  delegation  also  accepted 
French  invitations  to  visit  the  "Alian"  delivery  vehicle 
manufacturing  plant,  the  Paris  Aerospace  Exhibit,  the  French 
Aerospace  Research  Institute  in  Paris,  and  two  aviation  schools 
in  Toulouse. 


A  SUCCESSFUL  CONFERENCE  ON  TRANSONIC  AERODYNAMICS 

The  China  Aeodynamics  Research  Committee,  China  Aviation 
Committee  and  China  Aerospace  Institute  jointly  sponsored  a 
Transonic  Aerodynamics  Conference  in  Xian  from  May  26  to  30. 

Transonic  aerodynamics  includes  both  subsonic  and  supersonic 
aerodynamics  and  the  complex  phenomena  of  their  existing  shock 
boundary  layers  are  difficult  to  study  and  calculate.  Because  of 
the  application  of  electronic  computers,  research  work  on 
transonics  has  developed  quickly.  In  recent  years,  China’s 
aerodynamics  researchers  have  given  full  attention  to  the 
problems  of  this  area  and  have  attained  pleasing  results.  There 
were  74  professors,  specialists  and  engineers  from  31  units 
participating  in  this  academic  conference  and  42  papers  were 
read,  each  of  which  discussed  in  detail  technical  problems  in 
transonic  aerodynamics.  This  conference  had  the  effect  of  further 
promoting  Chinese  research  in  this  area. 
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THE  CSAA  SPONSORED  CONFERENCE  ON  GEARS  AND  REDUCERS  IN  AIRCRAFT 


The  CSAA  sponsored  a  conference  on  gears  and  reducers  in 
aircraft  in  Changsha  from  June  12  to  17  of  1981.  This  was  the 
first  national  special  conference  on  gears  and  reducers  in 
aircraft.  The  conference  received  65  papers  on  design,  technology, 
testing  and  measuring,  45  of  which  were  delivered.  Other  related 
specialists  were  invited  to  the  conference  to  give  special 
reports  on  the  developments  of  gears  and  reducers  in  aviation 
both  domestically  and  abroad. 

The  papers  received  by  the  conference  revealed  the  following 
outstanding  features:  firstly,  because  of  the  high  speed  and 
heavy  load  features  of  aviation  gears,  the  finite  element  method 
has  been  widely  applied  for  stress  analysis  so  as  to  raise  the 
gear's  strength  design  level.  Secondly,  research  on  the  design 
and  technology  of  the  height  and  shape  of  aviation  gear  teeth 
have  attained  new  advancements  as  well  as  certain  achievements 
in  reducing  gear  transmission  noise  and  vibration.  Furthermore, 
the  most  recent  research  achievements  on  the  development  and 
application  of  gear  integral  error  analysis  and  the  integral 
error  measurer  were  also  reported  at  the  conference.  These  topics 
drew  enthusiastic  interest  from  the  participants.  Several  papers 
systematically  introduced  the  newest  technology  recently  imported 
from  abroad  for  the  design,  manufacture,  testing  and  measuring 
of  gears  and  reducers  in  power  sets.  This  was  given  special 
attention  by  related  specialists  who  considered  many  of  the  new 
techniques  quite  worthy  of  study  and  research. 
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During  the  conference,  they  also  heard  an  introduction  on 
formulating  "criterion  for  gear  strength  in  aviation"  and  also 
held  discussions  on  this  special  topic.  All  of  the  participants 
considered  this  item  of  work  of  great  significance,  that  is, 
to  overcome  many  difficulties  to  attain  achievements  and  to  gain 
results  as  quickly  as  possible.  They  also  offered  valuable  and 
constructive  suggestions  for  future  work. 

To  sum  up,  the  participants  considered  that  this  conference 
reflected  the  new  levels  of  scientific  research,  design  and 
production  of  aviation  gears  and  reducers  in  China.  The  academic 
atmosphere  of  the  conference  was  enthusiastic,  the  representa¬ 
tives  spoke  out  freely  and  discussions  were  lively.  Anticipated 
goals  were  reached  and  excellent  goals  were  attained.  Finally, 
proposals  were  made  on  the  questions  of  the  arrangements  of 
future  activities  and  the  establishment  of  specialized  academic 
organizations . 
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FIRST  CONFERENCE  ON  SCIENTIFIC  MANAGEMENT  SPONSORED  BY  THE  CSAA 

The  CSAA  convened  the  First  Conference  on  Scientific 
Management  from  July  18  to  22  of  1981  in  Chengdu,  The  conference 
had  two  aims:  one  was  academic  exchange  and  the  other  was 
deliberations  on  the  setting  up  of  a  corresponding  academic 
organiztion.  Over  90  papers  were  collected  at  the  conference 
and  among  these  some  were  delivered  and  some  exchanged  in  small 
groups.  The  contents  of  the  papers  widely  reflected  the  basic 
research  developments  of  scientific  management  in  the  aviation 
industry  and  also  discussed  the  various  problems  of  organization, 
using  qualified  personnel,  economic  quality  and  systems  analysis, 
planned  appraisals  and  behavioral  science. 

The  conference  opened  with  enthusiasm  and  the  representatives 
spoke  freely  thus  realizing  an  academic  atmosphere  of  letting 
a  hundred  flowers  bloom  and  a  hundred  schools  contend.  Specially 
invited  representatives  gave  speeches  on  "Scientific  Policy- 
Making  and  Modern  Brain  Trusts,"  "Modern  Leadership  Techniques," 
:The  Application  of  the  Network  Method  in  the  Launching  of 
Chinese  Carrier  Rockets"  etc.  which  were  warmly  received  by  the 
representatives.  The  papers  acted  to  expand  everyone's  scope  and 
raise  their  understanding.  The  participants  further  comprehended 
that  scientific  policy-making  is  the  key  to  the  realization  of 
making  management  scientific,  is  an  important  factor  in 
guaranteeing  successful  development  in  social,  economic,  scientific 
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and  educational  areas,  and  is  also  a  basic  indication  of  testing 
the  level  of  leadership. 

The  conference  addressed  itself  to  the  problems  of  the 
development,  planned  management,  scientific  research  management 
and  economic  quality  of  China's  new  organizations,  and  held 
discussions  in  small  groups.  The  representatives  integrated  past 
model  developments,  introduced  special  research  experiences  and 
summarized  lessons.  It  was  unanimously  considered  that  model 
development  and  special  research  should  handle  well  scientific 
policy-making  and  also  should  use  the  network  method  to  guarantee 
coordination  between  the  various  systems.  In  this  way,  management 
will  become  scientific  and  raise  economic  achievements. 

To  sum  up,  after  participating  in  the  conference,  the  comrades 
understood  that  the  scientific  management  discussed  was  not  the 
general  planning  and  management  work  of  the  past  but  a  new 
branch  of  learning  which  integrates  the  social  sciences  and 
natural  sciences.  Therefore,  the  representatives  agreed  to  set 
up  a  Chinese  Aeronautics  Conference  Scientific  Management  Special 
Academic  Organization  to  further  promote  the  future  development 
of  management  science. 
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